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SENSITIVITY ANALYSIS AND OPTIMIZATION OF AERODNAMIC 
CONFIGURATIONS WITH BLEND SURFACES 


ABSTRACT 

' A novel (geometrical) parametrization procedure using solutions to a suitably 
chosen fourth order partial differential equation is used to define a class of airplane 
configurations. Inclusive in this definition are surface grids, volume grids, and grid 
sensitivity. The general airplane configuration has wing, fuselage, vertical tail and 
horizontal tail. The design variables are incorporated into the boundary conditions, and the 
solution is expressed as a Fourier series. The fuselage has circular cross section, and the 
radius is an algebraic function of four design parameters and an independent computational 
variable. Volume grids are obtained through an application of the Control Point Form 
method. A graphic interface software is developed which dynamically changes the surface 
of the airplane configuration with the change in input design variable. The software is 
made user friendly and is targeted towards the initial conceptual development of any 
aerodynamic configurations. Grid sensitivity with respect to surface design parameters and 
aerodynamic sensitivity coefficients based on potential flow is obtained using an Automatic 
Differentiation precompiler software tool ADIFOR. Aerodynamic shape optimization of the 
complete aircraft with twenty four design variables is performed. Unstructured and 
structured volume grids and Euler solutions are obtained with standard software to 
demonstrate the feasibility of the new surface definition. 
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Chapter 1 


INTRODUCTION 

1*1 Motivation 

Design and optimization of airplane components has become a primary ob- 
jective for most researchers in aerodynamic community. The sudden interest can be 
attributed to the introduction of complex and composite materials required by ad- 
vanced aerospace vehicles, such as National Aerospace Plane (NASP) and High Speed 
Civil Transport (HSCT) aircraft. Here, the interdisciplinary interactions are partic- 
ularly important because of extreme flight conditions. The design of such vehicles 
requires many analyses over a wide range of engineering disciplines. 

In the past, design of flight vehicles typically required the interaction of many 
technical disciplines over an extended period of time in a more or less sequential man- 
ner. At present, computer-automated discipline analyses and interactions offer the 
possibility of significantly shortening the design cycle time, while simultaneous mul- 
tidisciplinary design optimization (MDO) via formal sensitivity analysis (SA) holds 
the possibility of improved designs. Each analysis is based on solving mathematical 
models describing physical laws associated with a discipline. The mathematical mod- 
els are systems of algebraic, differential, or integral equations which are solved on 
discrete domains called "‘grids” on, around, and interior to the vehicle surface. The 
geometric requirements are the definition of the vehicle surface and the generation 
of grids onto which solutions of the mathematical models are obtained. In the opti- 
mization of aerospace-vehicle designs, engineering disciplines are interconnected and 
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affect one another. The effects can be realized in two ways: (1) The output from 
one discipline is the input to another. (2) The vehicle geometry changes in response 
to a discipline, therefore affecting other disciplines. In multidisciplinary analysis, the 
vehicle surface remains constant and all disciplines analyze their physics based on the 
same surface. Whereas, in multidisciplinary optimization the vehicle surface must be 
allowed to change. A complete design and optimization analysis using all the relevant 
disciplines is still a formidable task even for an isolated airplane component such as 
a wing or fuselage. The computational cost associated with such analysis can easily 
strain the capabilities of current supercomputers. The magnitude of this problem 
can be best appreciated when a discrete aerodynamic or structural design analysis 
can exhaust the computational capability of a medium size supercomputer. The un- 
derlying problem is the expensive cost of the analysis for each discipline involved. 
Clearly the aerodynamics involve non-linear physics and use of composite materials 
would require non-linear structural analysis as well. For a simple aeroelastic problem, 
the entire system matrix must be simultaneously solved using mostly implicit solvers. 
The extensive computational demand for such coupling of the governing equations, 
will likely limit MDO to only individual components such as a wing or wing-section. 
The cost of optimization operations are relatively small and manageable. Two gen- 
eral directions to overcome these difficulties have been proposed by different research 
groups. The first direction leads toward modifying the existing computational tools in 
order to obtain a relatively cheap and reliable technique for design and optimization. 
The usually favored direct solvers, with all their advantages, require extremely large 
computer storage even for 2D applications. 

Creating an airplane surface or any other object surface with design pa- 
rameters implies that there is an underlining set of rules or correspondences (model 
functions) that are driven by the parameters and independent computational vari- 
ables. Surfaces grids are discrete evaluations of the surface functions, and surface 



grids can be described as organized sets of points. Different discipline analyses and 
different techniques within a discipline most often require different grids to be gener- 
ated from the surface model. In an environment where the ability to quicklv change 
features of the geometry is nearly as important as the geometry itself, it is desirable: 
(1) to have the geometry model specified in terms of a small number of design pa- 
rameters; (2) to visualize the geometry and interact with it to explore the envelope 
of possibilities; and (3) to quickly extract grids and grid sensitivity for automated 
analysis (both low-level and high-level) and optimization. As the geometry becomes 
detailed, it is imperative that a CAD model, with its general characteristics be de- 
veloped, and any parameter-defined model should be upgraded with a conventional 
CAD system. Alternately, it would be desirable to incorporate a methodology like 
the one described here in a conventional CAD system. 

Design parameters can be classified according to whether or not they are cou- 
pled. Uncoupled design parameters influence the solution independently and would 
be the major contributors to optimization process. These parameters could be geo- 
metric, flow-dependent, or grid-dependent. The geometric design parameters specify 
the primary shape of a typical aerodynamic surface. Flow-dependent parameters are 
usually free-stream conditions such as free-stream Mach number or angle of attack. 
The grid-dependent parameters, relatively new in aerodynamic optimization, affect 
the interior and boundary grids; therefore, influencing the solution and optimiza- 
tion process. Traditionally, geometric parameters are considered the most affluent 
in aerodynamic optimization, although, optimization with respect to other design 
parameters is gaining respectability. For optimization with respect to geometric de- 
sign parameters, a perturbation in parameters affect the surface grid and the field 
grid which, in turn, affect the flow-field solution. There are two basic components in 
obtaining aerodynamic sensitivity. They are: (1) obtaining the sensitivity of the 



governing equations with respect to the state variables, and (2) obtaining the sensitiv- 
ity of the grid with respect to the design parameters. The sensitivity of the state vari- 
ables with respect to the design parameters are described by a set of linear-algebraic 
relations. These systems of equations can be solved directly by a LU decomposition 
of the coefficient matrix. This direct inversion procedure becomes extremely expen- 
sive as the problem dimension increases. A hybrid approach of an efficient banded 
matrix solver with influence of off-diagonal elements iterated can be implemented to 
overcome this difficulty. 


1.2 Literature Survey 

1.2.1 Aerodynamic Design and Surface Modelling 

Airplane design has historically been divided into three phases [l] 1 : concep- 
tual design, preliminary design, and detailed design. The conceptual design of an 
airplane usually begins with specifications for a proposed mission and rough sketches 
of the configuration. Geometry begins to evolve in the form of sets of connected 
points, and as the configuration approaches the end of the conceptual design phase, 
Computer-Aided Design (CAD) models are created. In the preliminary-design phase, 
high level analysis and testing of physical models are performed. Geometry for com- 
putational analysis and the construction of test models is extracted from the CAD 
model. In the detailed-design phase, the CAD model is the central design representa- 
tion, now containing detailed information for manufacturing the airplane. According 
to Raymer [2] design drawing is often carried out with a computer-aided drafting sys- 
tem where the aircraft geometry is represented by character-lines on its surface. This 
constitutes only a partial definition of the aircraft’s surface and the process of lofting 
between the character lines is required to create the complete aircraft surface. Thus 


'The numbers in brackets indicate references. 




there is a need for mathematical methods for representing or parametrizing curves 
and surfaces, which are flexible enough to represent a wide range of shapes in an easy 
and intuitive manner. It is also desirable to choose a method which uses few surface 
defining parameters so as not to overcomplicate the problem which would lead to 
an excessive use of computational time whilst at the same time to ensure sufficient 
flexibility in the surface in order to avoid trivial solutions [3]. 

One method of surface representation commonly used in computer-aided de- 
sign applications is that of Bezier surfaces [4], Here the defining parameters are the 
set of control points which form the characteristic polyhedron to which the surface 
then approximates. One advantage of this method is that the effect of changing a 
design parameter, i.e., the effect of moving a control point, on the surface shape is 
intuitively predictable. An improvement to this method is found in B-spline surfaces 
where each control point only influences the region of the surface close to it [5], Both 
these two properties are useful from the point of view of the end-user. 

By the late 1970s, the CAD/CAM industry recognized the need for a modeler 
that had a common internal method of representing and storing different geometric 
entities. At about the same time, three major groups looked at the possibility of using 
Non Uniform Rational B-Splines (NURBS). Boeing began developing the Tiger sys- 
tem in 1979. Integrating B-splines [6] with rational Bezier representations [7] quickly 
led to rational B-splines. SDRC (Structural Dynamics Research Corporation) pur- 
sued NURBS commercially and in 1978, the company started working on a modeler. 
The rapid proliferation of NURBS is due partly to their excellent properties and 
partly to their incorporation in such national and international standards as IGES 
[8], PHIGS+ [9], Product Data Exchange Specification, and International Standard 
Office, and Standard for the Exchange of Product Model Data. 

These methods, however, were not suitable to the problems investigated 
by Bloor and Wilson [10], since even the simple cubic Bezier surface had sixteen 
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control points each of which had three degrees of freedom. Also, the Bezier formu- 
lation was based on design by changing small regions of the surface independently 
whereas they were concerned with a more global approach to design. Bloor and Wil- 
son [11] introduced the method of generating free- form surfaces using solutions to a 
suitably choosen partial differential equation. By regarding a blend as a solution to 
a boundary-value problem and by choosing appropriate boundary conditions, they 
demonstrated that a solution to a suitably chosen elliptic PDE gave a smooth blend- 
ing surface that had the required degree of continuity with the primary surfaces to 
which it joined. Bloor and Wilson [12] have extended their work for approximating 
surfaces, which are the solutions of partial differential equations, in terms of B-splines 
so that they can be represented in a form compatible with more established surface 
design techniques. 

1.2.2 Grid Generation and Solution Methods 

In recent times techniques for the automatic generation of computational 
meshes have received much attention. This is primarily due to the fact that there 
has been an increased effort in the development of algorithms for the solution of the 
flowfield equations. Historically, many of the fundamental developements in theoreti- 
cal fluid dynamics have rested upon conformal mapping techniques for incompressible 
potential flow in which solutions on the boundaries can be obtained without resort to 
information in the field. Also panel methods [13], which utilize distribution of sources 
and siks on boundary surfaces, have played and continue to play an important role in 
aerodynamics. Recently, however, attention has been primarily focused on solution 
techniques for the Full Potential, Euler and Reynolds- Averaged Navier-Stokes equa- 
tions. These equations are formulated on the basis of the continum hypothesis. With 
computers restricted in memory and speed it is not possible to consider all points in 
the continum domain and hence it is necessary to select a subset of points within a 



domain at which flow quantities can be calculated. The combination of points and 
connections between points defines a mesh or grid on which numerical methods for 
the solution of the flow equations can be constructed. The assumption is then made 
that the information at these points is sufficient to describe the complete flowfield. 

In the most widely used approach [14] the domain is divided into a struc- 
tured assembly of quadrilateral cells. The structure in the mesh is apparent from the 
fact that each interior nodal point is surrounded by exactly the same number of mesh 
cells. Mesh generation, however, has proved to be a stubbornly difficult problem. 
Considerable effort has been devoted to this area in recent years as evidenced by the 
extensive literature on mesh generation. Numerical mesh generation techniques [15] 
have proved to be a powerful approach for creating meshes around complex shapes. 
Algebraic methods based on surface fitting [16], transfinite interpolation [17], and se- 
quential mapping [18] have also been applied to treat a variety of geometric shapes in 
both two and three dimensions. All of these methods, however, encounter difficulties 
when applied to complete aircraft configurations consisting of a wing, fuselage, tail 
and nacelles. A promising technique to tackle complex configurations is the use of 
a multiblock structure or a splitting-up of the space around the configuration into a 
number of smaller and topologically simpler regions. Separate meshes can be gen- 
erated for each block. In some cases [19], the mesh is required to blend smoothly 
together at block interfaces to provide a mesh that can be viewed as a single block 
by the flow solver. In other cases [20], the mesh is not required to connect smoothly 
at the interfaces and interpolation is needed to transfer flow information between 
separate blocks. Smith et al. [21] have generated grids around very complex config- 
urations and very promising results have been obtained. 

Neverthless, the generation of a mesh around a complete aircraft config- 
uration, including engine nacelles, has resisted the efforts of researchers until fairly 
recently. The first published calculations using a structured, conforming mesh around 
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a wing/ fuselage/nacelle/pylon combination is the work of Vigneron et al. [22] More 
recently, Sawada and Takanashi [23] generated a structured mesh to calculate the flow 
over a complete aircraft with wing mounted nacelles. These are striking successes in 
the generation of structured hexahedral meshes around complex configurations. 

The alternative approach is to divide the computational domain into an un- 
structured assembly of computational cells. The notable feature of an unstructured 
mesh is that the number of cells surrounding a typical interior node of the mesh is 
not necessarily constant. The nodes and the elements are numbered and, to get the 
information on the neighbours, we store the numbers of the nodes which belong to 
each element. There is no concept of directionality within a mesh of this type and 
that, therefore, solution techniques based upon this concept (e.g. ADI methods) will 
not be directly applicable. The methods which are normally adopted to generate un- 
structured triangular meshes are based upon either the Delauny [24] or the advancing 
front [25] approaches. Discretization methods for the equations of fluid flow which are 
based upon integral procedures, such as the finite volume or the finite element method, 
are natural candidates for use with unstructured meshes. The principal advantage 
of the unstructured approach is that it provides a very powerful tool for discretizing 
domains of complex shape [26,27]. In addition, unstructured mesh methods naturally 
offer the possibility of incorporating adaptivity [28]. Disadvantages which follow from 
adopting the unstructured grid approach are that the number of alternative solution 
algorithms is currently rather limited and that their computational implementation 
places large demands on both computer memory and CPU [29]. Further, these algo- 
rithms are rather sensitive to the quality of the grid which is being employed and so 
great care has to be taken in the generation process. 
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1.2.3 Sensitivity Analysis and Optimization 

Sensitivity analysis (SA) provides a natural systematic means for both an- 
alyzing and predicting the behavior of physical approximations and computational 
systems or for identifying significant input parameters in a system. The literature 
on sensitivity analysis and optimization is quite extensive. The pioneering work on 
sensitivity analysis for MDO was started from Sobieski [30] to the CFD community 
for extending their present capabilities to include sensitivity analysis of aerodynamic 
forces. Yates [31] developed an analytical approach using an implicit differentiation 
in combination with linearized lifting-surface theory to evaluate the sensitivity co- 
efficients. This can be used as a benchmark criteria for assessing the accuracy of 
approximate methods. Murthy and Kaza [32] developed a semi-analytical technique, 
using linear unsteady aerodynamics, to study an isolated wing-section and rotating 
propfan blades. Some aeroelastic analysis for transport wing has been investigated by 
Grossman et al. [33], where a coupled aerodynamic and structure model influences 
the design. Livine et al. [34] and a few other researchers focus on more complex 
interactions such as inclusions of active controls on the overall optimization process. 
A number of researchers have successfully pursued the quasianalytical approach to 
calculate sensitivity derivatives from nonlinear flow-analysis codes of varying degrees 
of complexity. For example, Elbana and Carlson [35] have computed wing-section 
aerodynamic sensitivity coefficients in transonic and supersonic flight regimes, and, 
more recently, they extended the work to 3D full potential equations using the sym- 
bolic manipulator MACSYMA to obtain the sensitivity coefficients. The procedure 
was applied to ONERA M6 wing platform with NACA 1406 wing sections [36]. The 
calculation of quasianalytical sensitivity derivatives is reported by Taylor et al. [37], 
Hou et al. [38], and Baysal et al. [39] for interior channel flows from a conven- 
tional upwind finite-volume solution strategy applied to the 2D Euler equations in 
body-oriented coordinates. These researchers have subsequently extended this work 
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to calculate sensitivity derivatives for 2D laminar flows from the thin-layer Navier- 
Stokes (TLNS) equations, including external flows over isolated airfoils [40]. Baysal 
and Eleshaky [41] presented an aerodynamic design strategy using direct differenti- 
ation of Euler equations. The procedure was applied to design a scramjet-afterbody 
configuration for an optimized axial thrust. This scheme was later extended to in- 
clude domain decomposition capabilities in order to reduce the computational costs 
associated with complex configurations [42]. Another strategy has been developed by 
Korivi et al. [43] and Newman et al. [44], where the sensitivity equations are recast 
and solved in incremental iterative form. The incremental iterative form is very flexi- 
ble and it increases the feasibility of solving the sensitivity equations for advanced 3D 
CFD codes. Korivi et al. [45] have demonstrated the use of this strategy to efficiently 
and accurately calculate quasianalytical sensitivity derivatives for a space-marching 
3D Euler code with supersonic flow over a blended wing-body configuration. 

Application of the quasianalytical methods requires the construction and 
evaluation of many derivatives, and for advanced CFD codes, the task of construct- 
ing exactly all of these required derivatives “by hand” is extremely complex. Ref- 
erence [43] shows that failure to consistently differentiate the turbulence modeling 
terms can result in unexpectedly large errors in the sensitivity derivatives that are 
calculated. A promising possible solution to this problem may be found in the use 
of a technique known as automatic differentiation (AD). Automatic differentiation 
is a chain-rule-based technique for evaluating the derivatives of functions defined by 
computer programs with respect to their input variables and has been investigated 
since 1960. Progress towards a general-purpose AD tool has been made with the 
development of ADIFOR by a joint effort of Argonne National Laboratory and Rice 
University. ADIFOR differentiates programs written in Fortran 77; that is, given 
a Fortran procedure (or collection of procedures) that describe a “function” and 
an indication of which variables in parameter lists or common blocks correspond to 
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“independent” and “dependent” variables with respect to differentiation, ADIFOR 
produces Fortran 77 code that computes the derivatives of the dependent variables 
with respect to the independent ones. ADIFOR has recently been tested by Bischof 
et al. [46] and Green et al. [47] in applications to an advanced CFD flow-analysis 
code called TLNS3D [48]. In these studies, a high Reynolds number, turbulent, 3D 
transonic flow over the ONERA M6 wing was selected as the example problem. 

1.3 Objectives of Present Study 

After reviewing relevant literature, it is aparent that parametrization of air- 
craft geometry plays an important role in the design process. Despite the differences 
in various approaches towards aircraft design it is agreed upon to identify an early 
stage in the design process during which general questions considering the aircraft’s 
configuration be studied; when, in order to meet whatever requirements exist, various 
alternative design solutions must be considered. In the past, when considering the 
question of the physical properties of new design, designers have had to rely upon 
their own knowledge and experience, and, further along in the design process, model 
testing. However, the increasing sophistication of numerical methods and the increas- 
ing power of computer hardware have meant that the properties of new design can 
be analysed by computer long before any physical model is created. Furthermore, 
whereas the main use of numerical methods has been an alternative to model test- 
ing, there is an increasing trend towards their use in the design process as a tool for 
optimization. Development of an efficient and reliable surface definition, grid gener- 
ation, grid sensitivity and optimization for conceptual design of aerodynamic shapes 
appears essential. 

An important ingredient of grid sensitivity and surface optimization is the 
surface parameterization. The most general parameterization would be to specify 
every grid point on the surface as a design parameter. This, although convenient, 



L*2 


is unacceptable clue to high computational cost. It is essential to keep the number 
of parameters as low as possible to avoid a surge on computational expenses. An 
analytical parameterization, may alleviate that problem but it suffers from lack of 
generality. A compromise would be using spline functions such as a Bezier or Non- 
Uniform Rational B-Spline (NURBS) function to represent the surface [49]. In this 
manner, most aerodynamically inclined surfaces can be represented with few control 
(design) parameters. This method has its own disadvantages, like the definition of 
wing fuselage intersection. The method of generating blend surfaces was a key area 
which led to the investigation of free form surfaces. Generation of free form blend 
surfaces was investigated by Bloor et al. [50] . The surfaces which they generated were 
quite interesting and the applications ranged from telephone handset to hull of a 
ship. They used the solution of fourth order partial differential equation to generate 
blend surfaces. In this study the idea was explored and was used towards generating 
aerodynamic shapes. 

With the advance in computers much research have been directed towards 
the development of graphic interface, which could accurately represent the surface on 
a computer screen. Most of the available graphics software have the capability of dy- 
namic translation and rotation. It was realized after reviewing the literature that the 
need for a graphic interface which could help the designer view the dynamic change 
in surface with the change in design variable was extremly helpful. This would act as 
an additional tool in the initial conceptual development of surfaces. 

The second main objective of this study is to do a grid sensitivity and surface opti- 
mization. Unlike aerodynamic considerations, the grid sensitivity analysis has been 
used on structural design models for a number of years. In this context, grid sen- 
sitivity can be thought as perturbation of structural loads, such as displacement or 
natural frequency, with respect to finite element grid point locations [51]. Two basic 
approaches have been cited for grid sensitivity derivatives. The first approach, known 
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as implicit differentiation, is based on implicit differentiation of discretized finite el- 
ement system. The other, which is based on the variation of continuum equations, 
is known as variational or material derivative approach. The main objective here 
is to develop a fast and inexpensive method for grid sensitivity to be used on an 
automated aerodynamic optimization cycle. Among two major classes of grid gener- 
ation systems (Algebraic, Differential), algebraic grid generation systems are ideallv 
suited for achieving this objective. The explicit formulation, resulting in a fast and 
suitable grid, enables direct differentiation of grid coordinates with respect to design 
parameters [52,53]. The development of software packages like ADIFOR, which could 
compute the derivatives in a manner that could save the time and effort of analytical 
methods was extremely helpful. This study involves the application of this software 
to compute both the grid and flow sensitivity towards an optimization study. 

The organization of this study is as follows. The physical and geometric 
representations of a typical model are derived in Chap. 2. Chapter 3 discusses the 
graphical user interface. The grid generation algorithm for both structured and un- 
structured is described in Chap. 4. The method of solution is provided in Chap. 5. 
Chapter 6 discusses the theoretical formulation and aerodynamic sensitivity equation. 
The results are presented and discussed in Chap. 7. Finally, some concluding remarks 
are provided in Chap. 8. 



Chapter 2 


PHYSICAL MODEL 

2.1 Computer-Aided Geometric Design (CAGD) 

In the late 1950s hardware became available that allowed the machining of 
3D shapes out of blocks of wood or steel [54]. These shapes could then be used as 
stamps and dies for products such as the hood of a car. The bottleneck in this produc- 
tion method was soon found to be the lack of adequate software. In order to machine 
a shape using a computer, it became necessary to produce a computer-compatible 
description of that surface. The most promising description method was soon identi- 
fied to be in terms of parametric surfaces. The theory of parametric surfaces was well 
understood in differential geometry. Their potential for the representation of surfaces 
in a Computer-Aided Design (CAD) environment were not known. The exploration 
of the use of parametric curves and surfaces to represent objects in computational 
environment [55] can be viewed as the origin of Computer Aided Geometric Design 
(CAGD). 

Surfaces can be defined by implicit algebric equations or explicit parametric- 
algeabric equation [56]. Parametric equations have dominated CAGD because of their 
intrinsic simplicity for modelling complex objects. 

In the development of parametric curves and surfaces, two different ap- 
proaches have evolved [57]. They are referred to here as “interpolative” and “approx- 
imative”. In an interpolative representation, points and derivatives on the curve or 
surface are used to control the formula defining the curve or surface. Langrangian and 
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Hermite interpolation formulas are examples of this approach. In an approximative 
approach points not necessarily on the curve or surface control the formula defining 
the curve or surface. The Bezier and B-Spline representations are examples of this 
approach. 

In the design process using an interactive CAD system, the approximative 
approach is highly advantageous. After prescribing an initial set of control points, 
the designer can pick and drag points and simultaneously observe the change in the 
shape of the surface. 


2.1.1 Geometric Representation of Wing Section 

The most commonly used approximative representation is the Non-Uniform 
Rational B-Spline (NURBS) function. The NURBS provide a powerful geometric tool 
for representing both analytic shapes (conics, quadrics, surfaces of revolution, etc.) 
and free- form surfaces [58]. The relation for a NURBS curve is 
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where X(r) is the vector valued surface coordinate in the r-direction, D, are the con- 
trol points (forming a control polygon), cu, are weights, and N lyP (r ) are the p-th degree 
B-Spline basis function defined recursively as 
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The r, are the so-called knots forming a uniform knot vector 
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where the end knots a and 6 are repeated with multiplicity p + 1. The degree, p, 
number of knots, m + 1, and number of control points, n + 1, are related bv 


m = n + p + 1. (2.4) 

For most practical applications the knot vector is normalized and the basis function 
is defined on the interval (a = 0, b — 1). Equation (2.1) can be rewritten as 


XH=I»)D, fl„„(r) = -g 1 




i = 0, ...., n 


>=o 


(2.5) 


Yli= 0 N;, p (r)u,’; 

where Ri, p (r) are the Rational Basis Functions, satisfying the the following properties 
among many others found in [59] 


Ri,p( r ) — 1 ^i,p( r ) — (“'6) 

1=0 

Figure 2.1 shows a six control point definition of the cambered airfoil ob- 
tained by Eq. (2.5). The points at the leading and trailing edge are fixed. Two 
control points at 0% chord are used to affect the bluntness of the section. The effect 
of the movement of the control points to create another airfoil is shown in Fig. 2.2. 
Figure 2.3 shows the effect of increasing the weight of the middle control point. It is 
seen that the curve is pulled towards the control point. An arc length distribution of 
the unit line is used for the knot vector 

An interactive program based on Eqs. (2. 1-2.5) have been developed. The 
program is menu driven, where after prescribing an initial set of control points, the 
designer can pick and drag these points and simultaneously observe the change in 
shape of the curve. Figure 2.4 shows the snap shot view of the interactive program. 
The cursor is drawn as a cross hair and different options are available in the pull down 
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Fig. 2.1 Six control point wing section definition. 
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Fig. 2.3 Effect of increasing the weight of control point 
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Fig. 2.4 Snapshot of graphic interface of NURBS curve. 
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menu. Weights associated with each of the control point can also be changed. The 
distribution of points for the NURBS curve is set to arc length formulation which can 
be also changed by th user. The program has the capability to output the NURBS 
curve into a predefined file. 


2.2 Surface Design and Parametrization 

The description of parametric surfaces in Computer-aided design can be 
broadly classified into the categories of shape representation and shape design [ 60 ]. 
Shape design is typically accomplished in an interactive manner i.e., the designer 
starts with a sketch and refines it untill it meets the requirements. Characteristically 
to the representation approach, there already exists a prototype of the model and 
numeric information that describes it, for which the corresponding computer model 
is processed automatically within the appropriate tolerance. Shape modifications 
have been of interest in both CAD/CAM and graphics for atleast two decades. In 
CAD/CAM, shape serves such purposes as aesthetics (free form design) , smooth- 
ing (removing wiggles, bumps etc.), satisfying special design requirements such as 
generating hard points or hard lines and adjustment of geometry (eg. spline based 
variational geometry). In graphics, shape can be used to generate a large variety of 
shapes or to perform animation based on subtle modifications. In any case, shapes 
generated by a computer system are rarely immediately acceptable, and subsequent 
modifications are required. The available techniques for modifications depend on the 
underlying representation scheme. Using splines, the modifications can be accom- 
plished by, moving control points (Refinement can be made by either degree elevation 
or knot insertion), using special blending functions such as tensioned splines, and 
using rational polynomials with weight. 
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For an aerospace vehicle such as High-Speed Civil Transport (HSCT), the 
traditional approach to design is for aerodynamics and performance disciplines to ini- 
tially create the vehicle surface[61]. The process is to define the planform, wing, fuse- 
lage, engine nacelles, and major control surfaces with aero/performance independent- 
design parameters. For instance, the wing is specified by the planform description, 
wing section, dihedral angles, and twist angles. Several sections are required for a 
wing. Approximately 50-100 independent parameters are required to specify a rough 
vehicle surface [62]. Usually a sparse set of points on the component surfaces which 
can be thought of as a very coarse grid becomes the surface description for analyses. A 
refined definition of a vehicle surface is obtained by applying Computer-Aided Design 
(CAD) techniques to sparse definition. The input to the CAD system is the sparse 
definition. CAD is used to create a patch definition of each vehicle components, and 
add surfaces such as fillets and wingtips. 

A patch is represented mathematically as 

m n 

%») = EEW(>W') (2-7) 

1=0 k -0 

0 < u,v < 1. 

where u and v are parametric coordinates, h is a matrix of surface definition param- 
eter and //-"(u), and H%(v) are interpolation functions respectively in the u and v 
directions. 

For the case of bicubic surfaces (m and n=3) the matrix of defining param- 
eters is 

ar(0,0) z v (0,0) x v (0, 1) a:(0, 1) 


A. , 


r* u (o,o) x u u(o,o) x u u(o, i) £ U (o, i) ■ 



[x u (l,0) x u u(l,0) x u u(l, 1) x u u(l,l). 


*(1,0) *v(l,0) x„(l,l) x(l,l) 

The elements of h lt k are the corner points of the patch derivatives with respect to 
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the parameter variables at the corner points anti cross derivatives with respect to the 
parametric variables at the corner points. For a bicubic patch there are 48 defining 
parameters, and a refined vehicle surface may consists of several hundred patches. 

In this study two methods of representing the vehicle surface are considered. 
The first is the most general approximative representation i.e., Non-Uniform Rational 
B-Spline [63] (NURBS) and the second is a novel parametrization procedure which 
uses the solution to a suitably choosen fourth order Partial Differential Equation [64] 
(PDE) to represent the surface. 

The commercial environment in which the two parametrization procedures 
was investigated requires that it should satisfy the following: 

(1) provide flexibility to design geometry 

(2) give a set of tools the designer can invoke at any stage of the design process 

(3) work in a reliable, fast and accurate manner 

(4) operate such that any modifications should preserve the entire continuity of the 
geometry, and 

(5) provide analytical equation defining surface to perform design optimization. 


2.2.1 M-6 Wing NURBS Representation 

A NURBS surface [65] is the rational generalization of the tensor product 
nonrational B-Spline surface and is defined as 


S(u,v) 




(2.9) 


E " =0 EJU N ‘A u ) N j,i( v )uij 

where are the weights, form a control net, and N hP (u) and Nj, q (v) are the 
normalized B-Splines of degree p and q in the u and v directions. The knot vectors are 


u = 


0, • • • , 0, Up_|_ i , ...., u r _ p— i, 1, ■ • • , 1 


P+1 


P+1 


( 2 . 10 ) 
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V ^ 0, • • • , 0, l^+X, Vj — rj— 1, !,•••, I 

' — v — ' > — v — ' 

l 9+1 q + 1 

where r = n + p + 1 and s = m + q + 1. 


Introducing the piecewise rational basis functions: 


the surface Eq. 


(2.9) can be written as 


Nj, p (u)Nj' q (v)LOjj 
E/lo Nk,p{u)Ni, r ,(o)u; k j 


( 2 . 11 ) 


(2.12) 


S(u,u) = £ 5Z (2-13) 

k = 0 j = 0 

A NURBS surface has the property £"=0 £"L 0 N t , p (u)N Jitl (v) = 1 and reverts 
to a B-spline when all the weights are 1. A NURBS surface has the advantage of 
being able to represent free-form surfaces, and with the proper choice of weights, 
conic surfaces. 

The surface skinning technique [66] is used to obtain the NURBS surface. 
The task of skinning is to fit a surface through an ordered set of space curves, called 
as section curves. The positioning of section curves in the three-dimensional space is 
customarily done with respect to a spine curve, from which appropriate orientation 
vectors can be automatically computed. The skinned surface is obtained in three 
steps: 

1. All the cross-sectional curves are first made compatible. That is, all the 
curves should have the same degree and number of control points and be defined over 
the same knot vector. 

2. Next u values and a knot vector V is calculated for interpolation with 
degree-q NURBS curves. 

3. Using the above values, curves are interpolated through the control points 
calculated by Eq. (2.13). 

ONERA M6 wing is used to demonstrate the skinning technique. The points 
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Fig. 2.6 Shaded NURBS suface 
for ONERA M6 wing. 


Fig. 2.7 Coarse surface grid over 
NURBS surface. 
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generated from the CAD software is used in the skinning process and the control 
points generated are shown in Fig. 2.5. The shaded NURBS surface is shown in Fig. 

2.6. A coarse CFD grid is generated on the surface of the wing and is shown in Fig. 

2.7. 

2.2.2 PDE Method 

The PDE method generates a surface X in Euclidean 3-space, which is a 
function of two parameters, i.e., X = (x(u,v), y(u,v), z(u,v)). The surface is obtained 
by solving a partial differential equation (PDE), in parameter u,v space, subject to 
boundary condition on X and its normal derivative with respect to u and v. In gen- 
eral, the order of PDE determines the number of derivatives of the unknown function 
that must be specified in the boundary condition. If control over both shapes of 
the curves bounding the PDE surface patch and the directions and magnitude of the 
coordinate vectors X u and X v at the edge of the patch are required then atleast a 
fourth order PDE is needed to generate the surface. The PDE may be written as 


d 2 


+ a 2 


d 2 


X = 0 


(2.14) 


L du 2 ' dv 2 l 
where X = (x(u,v), y(u,v), z(u,v)). 

The appropriate boundary conditions for Eq. (2.14) are the value of X and 
its normal derivative around the edges of the domain in the (u,v) plane. Since the 
generating equation, Eq. (2.14), is an elliptic PDE, the solution becomes very sensi- 
tive to the choice of boundary conditions. The boundary conditions act as a powerful 
tool for surface manipulation by a designer and can be used as a design parameter 
in an optimization process. The boundary conditions on function X are choosen that 
the curves forming the edges of the surface patch have the desired shape. The direc- 
tion of the vector X u and X u are tangential to the isoparametric lines on the surface. 
Therefore by altering the values specified for X' u and X v along the boundaries, one 
can effect the direction in which the surface moves away from the edges of the patch. 
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The general solution of Eq. (2.14) can be written in the form 

o 

* = *00 + EM n (ti)co$(nu) + B n (u)sin(anv)) (2.15) 

n = l 

where the coefficient function /l n ( u) and B n (v) are of the form 

A n {u) = a ni e anu + a n2 ue anu + a n 3 e~ a ™ + a n 4 «e“ antt 

B n (v) — b n \ e nv + b n 2 vc nv + 6 n3 e nu + b n ±ve 

The quantities a n i, a n2 , a n3 , a n4 , 6 nll 6 n2 , 6 n 3 ? 6 n4 are vector valued constants that 

can be found for a particular solution by Fourier analysis of the condition imposed 
on the isoparametric lines bounding the patch. 

Consider now, the problem of creating simple blends between two circular 
cross sections. For an illustrative purpose, consider the blend between a cylinder and 
a plane. It is necessary to set up the problem as a boundary value problem in (u,v) 
space with boundary conditions specified along curves in the (u,v) plane that corre- 
sponds to closed curves in E 3 . One of the boundary curve is taken to be the plan 
outline of the circular cylinder. Another boundary curve which is the definition of the 
plane is taken to be u = 1 and again is given parametrically in terms of v. Knowing 
that seperable solutions to Eq. (2.15) are of the form sinusoidal function multiplied 
by exponential function, the choice of boundary condition must reflect this. Thus, 
for this example, the flat plane is considered to be at z = 0 , and the curve is defined by 

x p = r p cos(v ) 
y P = r p sin(v) 

Zp = 0 


( 2 . 16 ) 
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This is the boundary condition on X that is applied at u = 0. Similarly, at u = 1, 
the curve for the cylinder is defined by 

x c = r c cos(v) 

Vc = r c sin(v ) 

z p = h (2.17) 

Since the generating equation is a biharmonic like partial differential equation, it re- 
quires derivative boundary conditions which are given at the plane by 

x' p = s 2 cos(u) 
y' p = s 2 sin(v) 

4 = ° (2-18) 

and at the cylinder by 

x' c = 0 

y' c = 0 

z' p = Sl (2.19) 

Figure (2.8) shows the blend between the circular cylinder and the plane. 
The different constants which act as design parameters are the radius of the cylinder, 
the height of the cylinder and the slope of the cylinder. For the plane it is the radius 
of the circular plane and its slope. The effect of varying the slopes of the cylinder 
and plane on the blend is shown in Figs. 2.8(b-d). In Fig. 2.8(b), the radius of the 
cylinder is very large compared to the plane and a very large slope is choosen for the 
cylinder. It is seen that the grid lines near the cylinder is orthogonal. In Fig. 2.8(c), 
the radius of the cylinder is reduced and also the slope of the plane. In Fig. 2.8(d) a 
negative value of the slope is choosen for the plane. 
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a) Regula Blend. 



b) Large radius and high positive 
value of slope for cylinder. 


c) Small value of slope for 
the plane. 



d) Negative value of slope for the plane. 


Fig. 2.6 Blend between a circular cylinder and a plane with change in 
design parameters. 
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It is concluded from these figures that the slope plays an important role in determin- 
ing the blend between two cross-sectional curves. 


2.3 Generation of Complete Aircraft by PDE Method 

Consider an aircraft shape made up of five patches: a fuselage, an inner 
wing, an outter wing and vertical and horizontal tails. For simplicity the fuselage 
is defined algebraically. The characteristic lines which form the boundaries between 
adjacent surface patches are 

(1) the curve where the inner wing meets the fuselage 

(2) the curve where the inner and outer wing meet 

(3) the curve at the tip of the outer wing and 

(4) the curves for the horizontal and vertical tails 

Figure 2.9 shows the different patches and sections used to represent the HSCT type 
configuration. 

A methodology based on the above mentioned theory has been developed to 
define a class of airplane configurations. It directly evaluates the surface grid, volume 
grid, and grid sensitivity, and the main objective of the methodology is to provide a 
grid generation package for conceptual design that could be used in a wide spectrum 
of analyses (potential flow to Navier-Stokes). The methodology and associated soft- 
ware is developed by Smith et al. [67] and is called Rapid Airplane Parametric Input 
Design (RAPID). 

The fuselage definition is an algebraic function which creates two surfaces - 
one above the fuselage intersection and one below. The airplane is considered to be 
symmetric about the x-z plane at y = 0, and only one side of the airplane is computed. 
The fuselage cross section is circular and is generated as a Fourier series whose axis 
is parallel to the x-axis, where the y and z coordinates of points on the surface are 
related by 
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x = RfCv = r (O cos ( 7r C/2), z = r(().sm(ff(/2) 

i/ 2 + ~ 2 = r 2 


( 2 . 20 ) 


with 


r (£) = a o sin(0) + a x sin{%0) 

where 


( 2 . 21 ) 


0 — 7r((l. — G^)^ + 02 ) 

In the preceeding equation a 0 , and a t are constants and 6 is a parameter which lie 
in the range 0 < 6 < 180. The value of £ = 0 corresponds to the end point on 
the fuselage, and £ = 0 corresponds to a point along the curve seperating the upper 
and lower fuselage surfaces. The parameters for the fuselage are: R F , the fuselage 
length; a 0 anda u control for the fuselage radius; and a 2 , a parameter to control a 
finite radius at the end of the fuselage. The boundary curve separating the upper 
and lower fuselage surfaces is a combination of the fuselage intersection with the lift- 
ing components and cubic curves connecting the intersections. The fuselage center is 
optionally allowed to translate upward along a quadratic function from the trailing 
wing/fuselage intersection point to the end of the fuselage. This creates a “duck tail” 
characteristic in the fuselage which can simulate take off and landing. Figure 2.10 
shows the fuselage cross section with different constants. 

2.3.1 Dirichlet Boundary Condition for the PDE Solution 

The curve where the outer wing and inner wing meet is taken as a plane 
curve (z=constant) having the shape of a simple airfoil. The airfoil shape at the crank 
is given by the relation 


x = Csin( 7rv) + X t 
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Fig. 2.9 Different patches and curves for the HSCT type 
representation. 
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<J Ucam T" ^ t 

z = o,q + H[ ( 2 . 22 ) 

where A f , Y t translates the crank boundary in a xy plane, also 


M 

Vcam ~ J^(2Lx - X 2 ) 


x < L 


(1 — 2Z. + 2£.c — x 2 ) 

Vcam — tri T \2 X < L 

1 1 — L) 


Dt = — —{sin(‘2xv) + / 3 1 sm(4jru) + P 2 sin(Q xv)) 


Parameters C, T, L and M are chord, thickness, location of maximum camber and 
maximum camber respectively, P l andP 2 are Fourier constants. The definition of the 
section starts at the trailing point, proceeds beneath the camber curve, around the 
leading point and over the camber curve back to the trailing edge. Figure 2.11 shows 
the airfoil definition at the intersection of the outer and the inner wings. 

The second character line lies on the surface of the fuselage. It is given 
parametrically by the equations 


Vj = y-Ta P + Y d 

zj = a(X) 2 -(y.T ap + Y d ) 2 (2.23) 

where B is the wing-root chord length, X d , Y d translates the wing fuselage intersection 
and T ap scales the thickness at the wing fuselage intersection relative to the thickness 
at the crank. This character line is basically a curve on the fuselage, whose projection 
onto the vertical plane containing the fuselage axis is an airfoil shape similar to the 
airfoil definition given by Eq. (2.23), but scaled by a factor (B/C). 

The third character line lies at the tip of the outer wing. It is given para- 


metrically by the equations 
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%ti — X c + 


X tl .x 

C 


M . _ x “-y , y 

Dll — Q + *C 


~ti — (lQ + H\ + H2 


(2.24) 


where X ti is the chord length at the wing tip, X c , Y c translates the wing tip in the xy 
plane, and H\. anilH 2 are the span length of the inboard and outboard wings respec- 
tively. 


The outer wing is generated by solving Eq. (2.14) using the boundary con- 
ditions obtained from the character lines, Eq. (2.22-2.24). Similarly the inner wing 
is generated with character lines given by Eq. (2.22) and solving Eq. (2.23). 


2.3.2 Neumann Boundary Conditions for PDE Solution 

Since the governing PDE equation, Eq. (2.14), is a fourth-order equation, it 
requires boundary condition on the normal derivatives of X (u,v) in the (u,v) parame- 
ter plane, which in the present case means boundary conditions on X u . The criterion 
used to decide how the boundary values of the tangent vector X u is chosen is based 
on the fact that, if tangent continuity between the blend and the primary surface is 
required, then the direction of surface normal must be continuous across the blend 
trimline. Note that the magnitude of this vector determines the ‘speed’ with which 
the isoparametric lines move away from the boundaries of the blend. 

On character-line (1), which lies at the junction of the outer and inner wing, 
the derivative boundary conditions are as follows: 

X n S j . X 

2 

Vu = 0 


z u — —S 1 


( 2 . 25 ) 



where S\ is an adjustable design variable. 


On character line (2) given by Eq. (2.25), which lies on the fuselage, the 
derivative boundary conditions are as lollows: 


, v q dy 

x f(u) = b 2 . — .sinivv 
Ov 

y f {u) = -S 2 .^.simrv 

_ jf - (Yd. + y.T ap )(-S 2 ^sinv) 

=J[U ~ a(Xy - (S 2 %sinvy 


(2.26) 


On the character line (3) which lies at the wing tip, the derivative boundary 
conditions are as follows: 


xu(u) = S 3 (— ^r^) 

yti(u) = 0 

zti(u) = 0 (2.27) 

The quantities S\, S 2 , and S 3 are adjustable design parameters whose values may 
be changed to control the transition of surface between inboard and outboard wing 
components and from wing into fuselage respectively. Figure 2.12 shows the complete 
PDE surface of the HSCT type configuration. Horizontal and vertical tails are added 
in a similar fashion as the inboard wing. All the major surface defining parameters 
are shown in the Fig. 2.13. 





Chapter 3 


GRAPHIC INTERFACE 

3.1 Introduction 

Interactive computer graphics is based on the concept of working with a 
model described by information stored in the computer. For a simple application such 
as drafting, the model includes only the information required to generate a picture 
of the physical object, such as the lines in a drawing or a detailed three dimensional 
representation. The application areas of simulation or computer-aided design and 
analysis involve a more extensive model in which the graphical data are associated 
with additional facts or mathematical equations explaining nonvisual characteristics 
of the physical object [68]. An abstract entity such as a chemical process can be 
modeled by a graphic flow chart. The work with the computer model is of two types, 
the creation of the model through input by the user and the display of the resulting 
model by the computer. The subject of interactive input encompasses more than just 
the way information is transferred from the input equipment to the graphics appli- 
cation program. In fact, the lowest level of input functions handle this transfer. It is 
the higher level functions that produce a satisfactory man-machine dialogue. During 
the design process the operator frequently needs to delete a previously drawn object, 
move an object, or modify an object in some way. The input commands that the 
operator can use to accomplish these actions all involve interacting with the drawing 
that is already stored in the data base. This is a more complicated process than sim- 
ply adding new data because it is first necessary to indicate to the computer which 
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Constructing user interfaces for programs is normally a time consuming 
process. They are though very important to help the user work with the program 
in an easy and pleasant way. In the past few years a large number of packages have 
appeared that help build up graphical user interfaces (so-called GUI’s) in a simple 
way. Most of them though are difficult to use and/or expensive to buy and/or limited 
in their possibilities. 


3.2 Interface Using Forms Library 

The Forms Library [69] package which was developed at the department of computer 
science, Utrecht University, Netherlands is a package that is simple to use, powerful, 
graphically good looking and easily extendable. 

The main notion in the Forms Library is that of a form. A form is a window 
(normally without a border) on which different objects are placed. Such a form is 
displayed and the user can interact with the different objects on the form to indicate 
their wishes. Many different classes of objects exist like for example, buttons that 
the user can push with the mouse, sliders with which the user can indicate a partic- 
ular setting, input fields in which the user can scroll through large amount of text, 
etc. Whenever the user changes the state of a particular object on one of the forms 
displayed the application program is notified and can take action accordingly. 

The forms library consists of a large number of C-routines and is simple to 
use. Defining a form takes a few lines of code and interaction is fully handled by 
the library routines. First one or more forms are defined, by indicating what object 
should be placed on them and where. After the form has been defined it is displayed 
on the screen and control is given to a library call fl-do-forms(). This routine takes 
care of the interaction between the user and the form and returns as soon as some 
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change occurs in the status of the form due to some user action. In this case control 
is returned to the program (indicating the object changed) and the program can take 
action accordingly, after which control is returned to the fl-do-forms() routine. Mul- 
tiple forms can be handled simultaneously by the system and can be combined with 
windows of the application program. 

An interface based on the forms library called Parametric Representation 
of Input Surface Mechanism (PRISM) is developed. The main application program 
where the points are generated to represent a surface is obtained from the program 
RAPID (67). The application program RAPID is converted to G' language and is 
combined with PRISM. The design variables which act as boundary condition for the 
solution of PDE equation are represented as different buttons on the screen. Buttons 
are provided for rotation and translation of the object and also to read and write a 
particular surface. The user can activate the program by simply typing PRISM. The 
program generates the surface points based on the application program which in this 
case is the solution of PDE equation. Each of the different sections is considered as a 
different surface and hence represented by a different color. As an example, a HSCT 
type configuration is considered and is represented by five different surfaces (fuselage, 
inner wing, outer wing, horizontal tail, and vertical tail). The user has the freedom to 
change the color of the different surfaces by selecting each of the surfaces individually. 
The program PRISM initially represents the surface in wireframe format which can 
be changed or rendered as shaded. Once the surface is represented on the viewing or 
main window, the user can pick any of the buttons of the different constants(design 
variables) and change to view the surface being changed interactively. The program 
runs in real time and gives a better understanding of the role played by each of the 
different design variables. A separate window is provided in PRISM which displays 
the numeric value of each of the design variables and also interactively shows the 
number being changed. The surface can be rotated and zoomed in and out. Once 
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the user is satisfied with the particular shape of the surface, it can be written out in 
a seperate file. Figure 3.1 shows the snapshot view of the software program PRISM. 
The wireframe surface mesh is shown with different values of the design variables 
in a seperate window. Figure 3.2 shows the change in airplane geometry when the 
parameters defining the wing, fuselage and grid concentration are changed. For con- 
vinence the program is also menu driven where the different options can be displayed 
by clicking the right mouse button. The rendering of surfaces in PRISM is done to 
better understand the curvature and roughness of the surface and is explained in Sec. 
3.3. 


3.3 Shaded-Image Rendering 

% 

High-resolution shaded raster images provide concrete visualizations of computer- 
generated surfaces. When features such as shadowing, specular reflection, and depth 
cueing are included, and the user is free to manipulate the viewpoint and the position- 
ing and intensity of the light sources, such images are extremly useful in understanding 
curved surfaces. 

The basic problem in generating high- resolution raster images is computing 
the intersections of a set of rays from the viewer’s eye with the surface. The approach 
to this problem depends on the surface formulation in use and the level of accuracy 
desired. 

A three-dimensional ray can be regarded as the intersection of two planes 

alx + 61 y + clz + dl = 0a2x + b2y + c2z + d2 = 0 (3.1) 

assumed to be nonparallel. The planes cut the surface in algebraic curves with poly- 
nomial equations of the form 
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Fig. 3.1 Snapshot of the interactive software PRISM. 
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Fig. 3.2 Snapshot of the software with the surface being modified. 
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Fl(u,v) = 0 , F2(u,v) = 0 

in the parameter space of the patch. These are obtained by direct substitution of the 
parametric surface equation r = r(u,v) into the plane equations given by Eq. (3.1). 
A simple lighting model for shading surface images typically includes ambient and 
directional light sources. The ambient component produces a uniform level of sur- 
face illumination, independent of viewing direction, while the directional components 
produce both diffuse and specular reflections with intensities depending on the angles 
between the surface normal and the viewing and illumination directions. 

Buttons are provided on the screen which when activated, render the surface 
of the geometry displayed as shaded. Colors are also provided for the user to choose 
for the surface rendering. A Toggle switch is provided which alternates the surface 
between rendering and wireframe. The surface is illuminated from a particular fixed 
direction which cannot be changed by the user. 


3.4 Interactive HSCT Shape Design 

To demonstrate the capabilities of the interactive program PRISM, a generic 
airplane shown in Fig. 3.3 is considered. The airplane defined with twentyone design 
variables and the values of each of these are shown in a separate window. The objec- 
tive of the transformation process is to obtain the HSCT type configuration shown 
in Fig. 3.4c, by interactively changing the values of the design variables. 

The design variables which are changed are ch; the wing chord length at 
mid section, XD; x-coordinate of wing trailing edge mid section YD; y-coorclinate of 
wing trailing edge mid section, HI; length of inner wing component, H2; length of 
outer wing component, B; chord length for root wing section, TAP; ratio of thickness 
at root to thickness at midsection, XTL; chord length of outboard wing section, XT; 
x-translation of outboard wing section, TL; length of fuselage, XTE; x-translation 
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of wing relative to fuselage, AO; maximum fuselage diameter, Al; fuselage tapper 
parameter and A2; parameter controlling end fuselage diameter. The snapshot of the 
different shapes attained during the interactive process is shown in Fig. 3.4. The 
method is extremely helpful in investigating the shapes for different airplane design. 

Four different airplane configurations generated with the help of PRISM are 
shown in Figs 3. 5-3. 8. Figure 3.5a shows the generic airplane defined by twentyone 
design variables. The software program PRISM has the capability to add or remove 
z-buffering and the user has the flexibility to adjust the direction of light. This fea- 
ture is captured and shown in Fig. 3.5b. Configuration 2 consists of moving the 
wing below the symmetry plane, giving the look of a high lift configuration. Here the 
fuselage diameter is increased and Fig. 3.6a shows the snapshot view. The position of 
wing and the ducktail fuselage is shown in Fig 3.6b. Configuration 3 which represents 
the HSCT type configuration is shown in Fig. 3.7a. This configuration is used in this 
study for a detailed analysis and results are presented in latter sections. Figure 3.7b 
shows the view from below. Figure 3.8a shows the fourth configuration made up of a 
delta wing, and in Fig. 3.8b a different view is shown. These figures not only shows 
the capability of the software PRISM, but also the flexibility of generating different 
airplane configurations using the PDE methodology. 
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Fig. 3.5a Generic airplane. 



Fig. 3.5b Generic airplane with different z-buffering 
and different intensity of light. 
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Fig. 3.6a Configuration 2 (low wing). Fig. 3.6b Configuration 2 showing 

the low wing. 



Fig. 3.7a Configuration 3(HSCT type). Fig. 3.7b Configuration 3 viewed 

from below. 







Chapter 4 


GRID GENERATION 

4.1 Introduction 

In recent times, techniques for the automatic generation of computational 
meshes have received much attention. This is primarily due to the fact that there has 
been an increased effort in the development of algorithms for the solution of the flow- 
field equations. Historically, many of the fundamental developments in the theoretical 
fluid dynamics have rested upon conformal mapping techniques for incompressible po- 
tential flow in which solutions on the boundaries can be obtained without resort to 
information in the field. Also panel methods, which utilize distribution of sources 
and sinks on boundary surfaces, have played and continue to play an important role 
in aerodynamics. Recently, however, attention has been primarily focused on solu- 
tion techniques for the Full Potential, Euler and Reynolds- Averaged Navier-Stokes 
equations. These equations are formulated on the basis of the continum hypothesis. 
With computers restricted in speed it is not possible to consider all points within a 
domain at which flow quantities can be calculated. The combination of points and 
connections between points defines a mesh or grid on which numerical methods for 
the solution of the flow equations can be constructed. The assumption is then made 
that the information at these points is sufficient to describe the complete flowfield. 

In order to study the flow-field around any aerodynamic configuration, a sys- 
tem of nonlinear partial differential equations must be solved over a highly complex 
geometry [70]. The domain of interest should be discretized into a set of points where 
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an implied rule specifies the connectivity of the points. This discretization, known 
as grid generation, is constrained by underlying physics, surface geometry, and the 
topology of the region where the solution is desired [71-72]. A poorly constructed 
grid with respect to any of the above constraints, may fail to reveal critical aspects 
of the true solution. 

The discretization of the field requires some organization in order for the 
solution to be efficient. The logistic structure of the data such as grid spacing, the 
location of outer boundaries, and the orthogonality can influence the nature of the 
solution [73]. Furthermore, the discretization must conform to the boundaries of the 
region in such a way that boundary condition can be accurately represented [74], 
This organization can be provided by a curvilinear coordinate system where the need 
for alignment with the boundary is reflected in routine choice of Cartesian coordinate 
system for rectangular region, cylindrical coordinate for circular region, etc. This 
curvilinear coordinate system covers the field and has coordinate lines coincident 
with all boundaries. To minimize the number of grid points required for a desired 
accuracy, the grid spacing should be smooth, with concentration in regions of high 
solution gradients. These regions may be the result of geometry (large surface slopes 
or corners), compressibility (entropy and shock layers), and viscosity (boundary and 
shear layers). A complex flow may contain a variety of such regions of various length 
scales, and often of unknown location. 

Two primary categories for arbitrary coordinate generation have been iden- 
tified. These are algebraic systems and partial differential systems. The algebraic sys- 
tems are mainly composed of interpolative schemes such as Transfinite Interpolation 
[75], Multi-Surface Interpolation [76], and Two- Boundary Interpolation techniques 
[77]. The basic mathematical structure of these methods are based on interpolation 
of the field values from the boundary. For partial differential equation systems, a set 
of partial differential equations must be solved to obtain the field values. 



The differential methods may be elliptic, parabolic, or hyperbolic, depending on the 
boundary specification of the problem. Each of these grid generation systems has 
its own advantages and drawbacks depending on geometry and application of the 
problem. Algebraic generating systems offer speed and simplicity while providing an 
explicit control of the physical grid shape and grid spacing. However, they might pro- 
duce skewed grids for boundaries with strong curvature or slope discontinuity. Partial 
differential systems, although offer relatively smooth grids for most applications, are 
computer intensive, specially for three-dimensional cases. An alternative, a common 
practice in recent years, has been to originate the grid using an algebraic system and 
then smooth the field using a differential system. Such hybrid approach has proven 
to be successful and cost effective for most applications. 

For complex geometries the multiblock mesh generation strategy is utilized. 
The idea behind multiblock mesh generation is that, instead of utilizing one global 
curvilinear coordinate system, several local curvilinear systems are constructed and 
connected together. The domain is subdivided into blocks and within each block a 
curvilinear system is derived. The block subdivision provides the necessary flexibility 
to construct structured meshes for complex geometrical shapes. The approach repre- 
sents a compromise between a globally structured mesh and an unstructured mesh. 

An array of general purpose grid generation softwares have emerged over the 
past few years. Among many others, the GRAPE2D of Sorenson [78], the EAGLE of 
Thompson [79], and GRIDGEN by Steinbrenner et al. [80] are the most widely used. 
The GRIDGEN series has both algebraic and differential generation capabilities on 
an interactive environment. The GRAPE2D solves the Poisson’s equation in two- 
dimension and utilizes a novel approach for determination of the boundary control 
functions. The EAGLE code combines techniques in surface grid generation as well 
as two or three-dimensional field grid generation. The ICEM/CFD has the capabil- 
ity of combining a full Computer Aided Design system (CAD), with grid generation 



module [81]. This provides an efficient and also quick procedure to reflect the CAD 
model changes on grids. Most of these packages furnish a host of options with a high 
degree of flexibility. However, intelligent use of the majority of these options requires 
the user to be well versed in current grid generation techniques. 

Over t he past few years, an alternative technique, unstructured tetrahe- 
dral grids, has received considerable attention [82]. In an unstructured mesh, unlike 
structured mesh, neighbouring points in the mesh in the physical space are not the 
neighbouring elements in the mesh point matrix. For any particular point, the con- 
nection with other points must be defined explicitly in the connectivity matrix. A 
constant reference to this matrix is made during the flow solution computation. In 
addition to their inherent capability of discretizing complex domain with ease, un- 
structured grids are suitable for efficient adaptive refinement, incorporation of moving 
boundaries, and local remeshing. These grids also offer better control over the mesh 
size and point distribution. In other words, unstructured grids are more flexible 
than their structured counterparts simply because of their irregularities. While in 
structured grids, mesh lines and planes should be continuous and conform to the 
boundaries and adjacent lines and planes throughout a domain, no such restriction 
exists in unstructured grids due to their lack of directionality. Generally, since tri- 
angles and tetrahedra are the simplest geometrical shapes having areas and volumes, 
respectively, they can discretize an irregularly shaped domain easier than quadrilat- 
erals and hexahedra. Furthermore, the number of neighbouring points surrounding 
each node in a structured grid is fixed, whereas in an unstructured mesh, this number 
varies from point to point. A consequence of this property of unstructured grids is 
that a large number of grid points on the surface of a geometry, where a fine resolution 
is required, do not have to be carried all the way to the outer boundaries where fewer 
points are needed. 
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There is a variety of methods for generation of unstructured grids in the liter- 
ature. Among these are Watson’s algorithm for Vornoi tessellation [83], the modified 
octree method [84] and the advancing front technique [85]. In this study, the software 
VGRID3D, a program for generation of three dimensional unstructured tetrahedral 
inviscid grids using the advancing front method [86] is used. This method is advo- 
cated here because it does not require a seperate library of modules to distribute 
grid points throughout the domain in advance like the Voronoi/Delauny family of 
unstructured grid generation techniques. 

4.2 Structured Grid Generation 

The majority of problems in physics and engineering can be described in 
terms of partial differential equations [87]. Many of these problems fall naturally into 
one of the three physical categories: equilibrium problems, eigenvalue problems and 
propogation problems. However, before solving such problems by numerical meth- 
ods, a system of partial differential equations should be solved to determine the mesh. 
The properties of meshes generated by this approach are intimately connected to the 
properties of the partial differential equations used as the mesh generation equations. 

Equilibrium problems are problems of steady state in which the equilibrium 
configuration is determined by solving a differential equation subject to boundary 
conditions. Such problems are known as boundary value problems and the governing 
equations for equilibrium problems are elliptic. 

Eigen value problems may be thought of as extensions of equilibrium prob- 
lems wherein critical values of certain parameters are to be determined in addition to 
the corresponding steady-state configurations. Propogation problems are initial value 
problems that have an unsteady state or transient nature. The problems involve the 
prediction of the subsequent behavior of a system given the initial state. The 
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governing equations for propogation problems are parabolic or hyperbolic. 

Structured algebraic grid generation techniques can be thought of as trans- 
formation from a rectangular computational domain to an arbitrarily shaped physical 
domain as shown in Fig. 4.1 [88]. The transformation is governed by vector of control 
parameters, P, and can be expressed as 


where 


xtf.if.c, P) 


'*K,*,c,pr 

< y(C»/>C»P) > 

U(c»7,c,p) J 


(4.1) 


0 < < 1, 0 < T] < 1, and 0 < C < 1. 


The control parameter P, is composed of parameters which control the primary 
shape of the boundary (design parameters), and parameters which control the grid 
(grid parameters). A discrete subset of the vector- valued function X(£, 77 j, P) = 
X{ X y z Y i jk = X* is a structured grid for & = = j^rj,Cfc = 7tEt> where 

i = 1, 2, 3 • • • , L, j = 1,2, 3, • • • , M and k = 1, 2, 3, ■ • • , N. 

Surface mesh generation is one of the most difficult and yet important as- 
pects of the total mesh generation problem. The surface mesh influences the field 
mesh close to the configuration, where flow gradients are important and need to be 
resolved accurately. Surface meshes have the same requirements for smoothness and 
continuity as the field meshes for which they act as boundary conditions, but, in 
addition, they are required to conform to the configuration surfaces, including, lines 
of component intersection, and to model regions of high surface curvature. 

In the software program RAPID geometric surfaces are generated using par- 
tial differential equation described in Sec. (2.6). The surface grid is created by 
evaluating the surface functions at discrete £(I) and £(K). In order to concentrate the 
grid in certain regions, such as wing/fuselage intersection, it is necessary to create 
control functions that map 0 < £,£ < 1 into 0 < £,£ < 1. The spacing of grid points 
within the topology constraints is very important for achieving acceptable accuracy 
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(a) Physical domain (b) Computational domain 


Fig. 4.1 Physical and computational coordinates. 
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in the application of a flow analysis about the vehicle surface. A double exponential 
function [89] which maps the computational variables £,77, and ( onto themselves is 
used here. The grid spacing control function is expressed as 




Figure 4.2 is used to help describe the grid control parameters AT, AT, A' 3) 
and A' 4 . Parameters AT and AT are coordinates of a point in the unit square. The 
quantity v is the independent computational variable and corresponds to the per- 
centage of grid points in a particular direction. The quantity v is the dependent 
computational variable and corresponds to the percentage of distance in the physical 
space along a grid curve. The parameters AT and A’4 are coefficients in the expo- 
nential functions defined for a particular part of the unit square. Where there is low 
slope in the control functions, there is a concentration in the grid points, and where 
there is high slope, there is dispersion in the grid points. In the RAPID methodology, 
Eq. (4.1) is used several times. The approach specifies a desired spacings at the 
v — 0 and/or at v = 1 and/or A'3. The parameters AT, A2, and AT are determined 
by a Newton- Raphson process while satisfying a first derivative continuity condition 
at (AT, AT). Figure 4.3 shows the grid distribution achieved on the fuselage by using 
Eq. (4.1) and Fig. 4.4 shows the grid distribution on the wing and wing fuselage 
intersection. 


The grid control parameters are distinguished from the configuration design 
parameters. The design parameters are referred to as the set V, and the grid 





Fig. 4.4 Grid distribution on the surface of the wing and the nose 
of the fuselage. 
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parameters are referred to as the set 1C. 1C which includes the grid spacing parameters 
described above and the volume grid control points discussed in the next section. 

4.2.1 Boundary Discretization and Volume Grid Around The Airplane Ge- 
ometry 

The orientation of the computational coordinates relative to physical coordi- 
nate, known as grid topology, is an important aspect of the transformation proceduie. 

In order to establish a grid topology for any geometry, it is essential to examine each 
component separately [90]. For any given geometry, there are several possible topolo- 
gies with different characteristics in terms of efficiency, coordinate cuts, singularities, 
etc. For example a typical wing-section geometry, may have at least three types of 
different topologies (e.g., C-, 0-, or H-types). The C- and O-type topologies usu- 
ally produce the most efficient grid. This topology produces no singularity and it is 
relatively simple to implement. For wing-sections with sharp noses, a H-type topol- 
ogy would be more appropriate. For more complex geometries, selection of different 
computational coordinate systems for different regions of physical domain might be 
required. In this case, physical domain is mapped into several computational sub- 
domains, where each sub-domain is reffered as a block. Therefore, it is possible to 
have a boundary-fitted coordinate system for a highly complex configurations. For 
the present study, the airplane geometry consists of two main components: the fuse- 
lage and the wing. The fuselage has a circular like cross-section which suggests that 
a natural O-type (cylindrical coordinate) grids. This topology produces a nearly 
orthogonal grid with one line polar singularity at the nose. For the streamwise direc- 
tion, it is feasible to have either a C-type or a H-type grid depending on the slope 
of the nose. For a fuselage with small nose slope, a H-type grid in the streamwise 
direction would be more appropriate. A wing has its own natural coordinates which 
are usually not compatible with the fuselage’s coordinate system. It is possible to 
generate a H-, 0-, or a C-type grids in the streamwise direction, and a C- or a H-type 
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in crosswise direction. To maintain a minimum of C° continuity at the interfaces, it 
is essential to select a compatible topology for the wing and fuselage. For most cases 
it is conceivable to generate a single block grid about these components, but this grid 
tends to be skewed for any practical purposes. A dual-block grid possesses much less 
skewness than a single-block grid. It consists of two large blocks, one covering the 
top portion of the physical domain, and the other covering the bottom portion of the 
physical domain. The dual-block topology is a direct consequence of using a H-type 
grid for the wing of zero wing-tip area. Figure 4.5 illustrates the mapping of a generic 
airplane geometry using a dual-block topology. A C-0 type grid have been chosen for 
a fuselage while the wing, horizontal, and vertical tails mapped to a H-H type grids. 

A Control Point Form/Transfinite Interpolation technique[91] is used to 
compute volume grids for the RAPID methodology. A considerable amount of in- 
formation has been published on this grid generation method and its variations, and 
only the major steps are presented here. 

Having established a grid on the configuration surface, the volume grid gen- 
eration is accomplished in four major steps described below. 

Step 1 is the determination of a grid in the symmetry plane. The basic 
functions used in RAPID are those for Bezier curves computed with the de Casteljau 
scheme[92]. Control points for an intermediate curve and for a far-field curve are com- 
puted from the dimensions of the fuselage, Fig. 4.6 . A set of points are distributed in 
the {-direction on the control curves obtained from the control points. Interpolation 
from the fuselage surface across the control curves is obtained with a de Casteljau 
application in the //-computational direction, and Fig. 4.7 shows the symmetry grid. 

Step 2 is the determination of a three-dimensional grid surface containing 
the lifting components shown in Fig. 4.8. Note that in the H-topology, the top 



62 


and bottom grids are considered separately. A process similar to that used with the 
symmetry grid for computing control points from the fuselage and lifting surfaces is 
applied. 


Step 3 is the determination of a cap grid. Control points are extracted from 
the extreme x and y grid coordinates in the lifting surface grid and the extreme z-grid 
coordinates in the symmetry plane grid. This is shown in Fig. 4.9. Casteljau scheme 
is applied with these control points, and the outer grid surface is shown in Fig. 4.10. 

Step 4 is the application of Transfinite Interpolation to compute the inte- 
lioi grid. Figure 4.11 shows a sample grid around the HSCT type configuration. . 

It is necessary to use several grid-spacing control functions and their control 
parameters in addition to the interpolation control points in order to achieve a good 
grid for a given set of design parameters. This requires some trial and error before ac- 
ceptable parameters are realized. However, once an acceptable set of grid parameters 
fC is found for a given set of design parameters V, small changes in V do not require 
changes in K. Therefore, repetitive small changes in the design parameters such as 
during configuration optimization, do not require the constant modification of the 
grid parameters. Also note that the volume grids obtained with this algorithm are 
computed only out to the wing tip. An additional far-field grid would be necessary 
for most high-level fluid analyses. 

A complete volume grid which extends beyond the tip of the wing surface is 
computed by GRIDGEN software. A comparitive study of the grids generated from 
RAPID is made with standard grid generation software GRIDGEN. Among the dif- 
ferent softwares available, the GRIDGEN software developed by MDA Engineering 
is used to develop grids around surfaces generated from RAPID. 




Fig 4.5 Dual-block grid topology for a generic airplane. 





Fig. 4.6 Symmetry plane control net. Fig. 4.7 Symmetry grid. 



Fig. 4.8 Surface grid containing lifting Components. 



Fig. 4.9 Control point for outer grid surface. 
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Fig. 4.10 Outer Grid Surface. 


67 



Fig. 4.11 Sample Grid Surfaces. 


\, v *i** 
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Fig. 4.12 Volume grid around the PDE surface using GRIDGEN. 
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The GRIDGEN software consists of three main modules namely, GRID- 
BLOCK, GRIDGEN2D, and GRIDGEN3D. The GRIDBLOCK begins with a dis- 
play of the 3D database of the airplane configuration. The 3D lines representing the 
bounding edges of the blocks are drawn. Once several connectors are added to the 
system, they are grouped together and assigned to blocks. Then computational direc- 
tions and dimensions on that block are defined. In the end flow boundary conditions 
and interblock connections are determined and assigned. The GRIDGEN2D is used 
to generate grid on the edges and surfaces. There are five modes of elliptic solvers in 
GRIDGEN2D. The first three solve directly for the Cartesian grid point coordinates 
in an iterative process and the next two solve the grid in parametric coordinates. 
This surface grid generation procedure is repeated for each surface in the face, and 
for each face in the block. The third and final step of the grid generation process is the 
distribution of grid points within the interior of each block. This task is performed 
with the batch code of GRIDGEN3D. Figure 4.12 shows the volume grid generated 
around the PDE surface. 


4.3 Unstructured Grids 

One of the greatest concerns in computational fluid dynamics is the gen- 
eration of suitable grids. Although considerable effort has been devoted towards 
development of robust and automatic grid generation methods, the process of gener- 
ating 3-D grids around complex geometries remains a formidable challenge. With the 
availability of large supercomputers, it is now possible to compute flowfields around 
complex configurations in a matter of hours. However, the process of grid genera- 
tion, using conventional structured grid methods, still makes up a large portion of a 
typical computational effort for a complex configuration. Use of unstructured grids 
has grown considerably in recent years due to their ability to produce quality grids 
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around complex configurations with relative ease. 

In recent years a wide variety of algorithms has been devised for the gen- 
eration of unstructured grids around bodies of complex geometrical shapes. Among 
the different techniques are the Watson’s algorithm for Vornoi tesselations, the mod- 
ified octree method [94] and the advancing front technique. Baker’s implementation 
and optimization of the Vornoi algorithm [93] has shown that fast and reliable grid 
generators for tetrahedral meshes can be produced. In this study, advancing front 
technique is used for grid generation, because it can easily be used for grid generation 
with directional refinement. Also it does not require a separate module to distribute 
points like the Delauny triangulation. 

4.3.1 Advancing Front Technique 

In the advancing front method, a grid is generated starting from the domain 
boundaries marching towards the interior of the computational domain. Unlike De- 
launy triangulation technique in which grid points are first distributed in the entire 
field and then connected to form cells, an advancing front introduces new points to 
the domain as tetrahedrons are made. The configuration of interest is first defined in 
terms of a number of surface patches. These patches are then triangulated to form 
the initial front. The front is then projected to the original surface which in this 
case is the NURBS and PDE surface. Next, tetrahedral cells are generated on top 
of triangular faces on the front by introducing new or using existing points. During 
this process, old faces are replaced by new ones, and the front is advanced in the field 
until the whole region is filled with grid cells. 

The entire grid generation process is summarized in the following main steps: 

(a) The boundaries of the domain to be grided are divided into a number of surface 
patches. These surfaces define the configuration of interest as well as the far-field 
boundaries. 

(b) A background grid is set up to define the local grid characteristics such as grid 
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point spacing. The spacing interpolation in the VGRID system is based on a struc- 
tured background grid [94]. This technique simplifies the specifications of grid density 
by introducing nodal and linear sources. The contribution from nodal sources are in- 
versely proportional to the square of the distance and the contribution from the linear 
sources are modeled similar to the diffusion equation. 

(c) Each surface patch is, in turn, subdivided into a number of triangles to form the 
first front (surface grid). 

(d) The triangles are then projected on to the actual surface which in the present 
case is the NURBS and the PDE surface. 

(e) The front is advanced in the field by introducing new points and forming tetra- 
hedras and new faces to complete the grid. 

(f) The completed grid may optionally be post-processed. 

The above described procedure is applied to an ONERA M6 wing. The wing 
has a leading edge sweep of 30 degrees, an aspect ratio of 3.8, a taper ratio of 0.56, 
and symmetrical airfoil sections. The wing has a root chord of 0.67 and a semispan 
b of 1.0 with a rounded tip. The computational domain is bounded by a rectangular 
box with boundaries at 

-6.5 < x < 11.0,0.0 < y < 2. Sand — 6.5 < * < 6.5 

Figure 4.13 shows the ONERA M6 wing bounded by the rectangular box. The M6 
wing is attached to one of the surface of the box. Triangulations starts from the 
surface of the box and the M6 surface and proceeds towards the interior of the domain. 
Figure 4.14 shows the surface triangulation on the actual NURBS M6 surface. This 
triangulation is obtained by projecting the initial triangulation of the surface on to 
the actual NURBS surface. 

Surfaces obtained from RAPID is also triangulated. In this case the HSCT 
type configuration is placed in the middle of the rectangular box. Figure 4.15 shows 
the surface mesh with the rectangular box and the HSCT type configuration. Figure 
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4.16 shows the surface mesh without the horizontal and vertical tails. In order to 
simulate the configuration with engines, two tapered rectangular boxes are placed 
just below the wings, and Fig. 4.17 shows the surface mesh. Figure 4.18 shows the 
surface mesh with horizontal and vertical tails. 
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Fig. 4.13 Far Field boundary for the ONERA M6 wing. 
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Fig. 4.14 Surface mesh on the ONERA M6 wing. 
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Fig. 4.15 Far field boundary for the HSCT type configuration. 




Chapter 5 


GOVERNING EQUATIONS FOR FLOW 

SOLUTIONS 


5.1 Unstrurctured Grid Solution 

The inviscid flow field is computed on the unstructured grids using USM3D, 
a three-dimensional upwind flow solver developed at NASA/LaRC [95]. The fluid 
motion is governed by the time dependent Euler equations for an ideal gas which 
express the conservation of mass, momentum, and energy for a compressible inviscid 
nonconducing fluid in the absence of external forces. The equations are given below 
in integral form, for a bounded domain fl with a boundary dQ, are expressed as 

l/// n Qdv+ /L F ( Q)fids=0 < 51 > 

where 

P 

pu 

Q = < pv > 

pw 

. e 0 . 

and 

p o 

pu ri x 

F(Q) • n = (V • n) < pv > + p < ri y > 

pw n 2 

, e 0 + P J o ; 

The equations are nondimensionalized with reference density p ^ and a speed of sound 
cioo. Here n x , n y and ri z are the Cartesian components of the exterior surface unit 
normal n on the boundary dCl. The Cartesian velocity components are u, v, and 


78 



79 


w in the x, y and z directions, respectively. The term e 0 is the total energy per 
unit volume. With the ideal gas assumption, the pressure and total enthalpy can be 
expressed as 

P = (7 - 1 ) (e 0 - i/> (u 2 + v 2 + tu 2 )^ 

h 0 = — ^ (u 2 + v 2 + w 2 ) 

7 — 1 p 2 v y 

where 7 is the ratio of specific heats and is prescribed as 1.4 for air. 

The spatial discretization is accomplished with a cell centered finite volume 
formulation using the flux difference splitting procedure. The solution is advanced 
in time using a three-stage Runge-Kutta time stepping scheme. Local time stepping 
and implicit residual smoothing are used to accelerate the convergence of the solution 
to a steady state. 

Boundary Conditions 

For the solid boundaries such as the wing and centerplane, the flow tan- 
gency condition is imposed by setting the velocities on the boundary faces to their 
cell center values and then subtracting the component normal to the solid surface. 
Density and pressure boundary conditions are simply set to the cell-centered value. A 
condition of zero mass and energy flux through the surface is ensured by setting the 
left and right states of solid boundary faces equal to the boundary conditions prior 
to computing the fluxes with Roe’s approximate Riemann solver. 

Characteristic boundary conditions are applied to the far-field subsonic boundary 
using the fixed and extrapolated Riemann invariants corresponding to the incom- 
ing and outgoing waves. The incoming Riemann invariant is determined from the 
freestream flow and the outgoing invariant is extrapolated from the interior domain. 
At an outflow boundary, the two tangential velocity components and the entropy are 
extrapolated from the interior, while at the inflow boundary they are specified as 
having far-field values. 
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The unstructured grid generated around the NURBS M6 wing using the 
advancing front technique (shown in Fig. 4.13) is used as a test case. The solutions 
were started from freestream initial conditions with first-order scheme until the L 2 - 
norm (RMS average of all residuals) decreased one order of magnitude, at which time 
the solver automatically switched to a higher order scheme. Converged solution is 
obtained for A/,*, = 0.84 and a — 3.06°. The upper surface pressure contour with 
contour intervals of A ( — 2 — ) = 0.02 is shown in Fig. 5.1. The figure clearly shows 
a double shock wave on the upper surface and is in good agreement with the results 
obtained by Frink et al. [96]. 

Converged solutions are also obtained for HSC'T type PDE surface shown in 
Figs. 4.16-18. Converged solutions are obtained for A/ t30 = 0.84 and a = 5°. Figure 
5.2 shows the shaded Cp plot. Contours are plotted by taking a cutting plane at the 
mid section of the configuration. A shock wave is seen at the upper surface of the 
wing. A total lift of 0.33358 and a drag of 0.04301 were obtained. To simulate this 
HSCT configuration with engines and to study the performance features, two engines 
of tappered square cross-section are placed below the wing. The unstructured grid 
shown in Fig. 4.17 is used for this case. Figure 5.3 shows the shaded Cp plot and for 
this case. A total lift of 0.313434 and a drag of 0.05932 were obtained. 

5.2 Potential Flow Solution 

A low-order potential-flow panel code for modeling complex three-dimensional 
geometries is used to calculate surface pressure variations. The flow field is assumed 
to be inviscid, irrotational and incompressible. The velocity potential is given by the 
Laplace’s equation: 

V 2 $ = 0. (5.2) 

The potential at any point P may be evaluated by Green’s Theorem which results in 
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p 0.812709 
g 0.718151 
S 0.625592 
0.529034 


Fig. 5.2 Cp plot for the HSCT configuration without engines. 




Fig. 5.3 Shaded Cp plot for the HSCT type configuration 
with engines. 
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the following integral equation 


*'-hJL 


S+IV+So 


(f - 1>,) -i ■ V (i) dS - ± jj s n (v« - V«.) «/S (5.3) 


It is assumed that the wake is thin and there is no entrainment, so the source term for 
the wake disappears and the jump in normal velocity across the wake is zero. Hence 


the simplified equation becomes 


* p = b //* 1 { * - v (?) dS - i JI S s ■ c v* - v*.) ds 

+ s JL { * u - * - l) " ■ v (?) (?) ,is + +■** < 5 - 4 > 

The Dirichlet type boundary condition is used to solve Eq. (5.3). The total 
potential $ can be viewed as being made up of an onset potential $oo and a pertur- 
bation potential <j> = $ — The potential of the fictious flow is set equal to the 
onset potential, <f> qq. With this boundary condition, the singularities on the surface 
tend to be smaller than if the potential of the fictious flow is set to zero because 
the singularities only have to provide the perturbation potential instead of the total 
potential. The general equation for the potential at any point P can be written as 

= l //s. p ' , "' v (?) «+**!+//, (?) dS +l V (?) (5.5) 

where I< = 0 if P is not on the surafec, I< = 2zr if P is on a smooth part of the 
outer surface, and K = -2? r if P is on a smooth part of the inner surface. If the sur- 
face is broken up into panels, Eq. (5.4) can be written in discretized form, breaking 
the integrals up into surface integrals over each panel. A constant strength source 
and doublet distribution is assumed over each panel and so the doublet and source 
strengths are factored out of the integrals. Taking point P to be at the centroid on 
the inside of one of the panels, the surface integrals over each panel are summed for 
all panels. For the panel containing point P, the surface integral is zero and only 



the — 27T/( p term remains in the bracketed part of Eq. (5.4). For all other panels, 
the surface integral is used and the —27 vf.i p term is zero since the pont P is not on 
the surface of any other panels. The process is repeated for point P at the centroid 
of every panel to yield a set of linear simultaneous equations to be solved for the 
unknown doublet strength on each panel. The surface integrals represent the velocity 
potential influence coefficients per unit singularity strength for panel K acting on the 
control point of panel J. Hence Eq. (5.4) becomes 


Ns 


N s 


N w 


Y (pkCjk) + Y ( a i<Bji<) + Y (pwlCjl) = 0 

K = 1 A'=l K = 1 


(5.6) 


where 


Bjk = If l -dS (5.7) 

J JK V 

and 

Cj K = [[ n-V-JS (5.8) 

J JK r 

The coefficients Cjk and Bjk represent the velocity potential influence coefficients 
per unit singularity strength for panel K acting on the control point of panel J. 
Equations (5.6) and (5.7) are functions of geometry only and thus can be solved for 
all panels to form the influence coefficient matrix. Since the source values are known, 
they may be transferred to the right hand side of the matrix equation. Solutions for 
Eqs. (5.6) and (5.7) can be found in [97]. 

As a test case, the PDE surface shown in Fig. 2.12 is considered. Only half 
of the configuration was modeled in PMARC. The other half of the configuration 
was generated by reflecting the model across the plane of symmetry. The wing was 
represented with 300 panels: 15 divisions in the chordwise direction on the upper and 
lower surfaces of the wing with denser spacing near the leading and trailing edges, 
and 10 divisions in the spanwise direction with denser spacing near the root and tip 
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of the wing. The tip of the wing was closed off with a flat tip patch. The fuselage 
was represented with 320 panels. The wing/fuselage junction was modeled such that 
wing and fuselage panels matched up exactly. An initial wake was attached to the 
trailing edge of the wing and to the aft of the fuselage and carried downstream 20 
chord length. Three time steps were specified to allow the wake start to roll up. The 
model was tested at an angle of attack of 4°, and the Cp plot is shown in Fig. 5.4. 



Fig. 5.4 Shaded Cp plot for the potential flow. 



Chapter 6 


METHOD OF SOLUTION FOR SENSITIVITY 

EQUATIONS 

6.1 Introduction 

Computational Fluid Dynamics (CFD) is now routinely applied to simulate 
flow about aerodynamic configurations. On current supercomputers, these simula- 
tions can require several hours per steady-state solution for viscous-compressible flow 
about airplane configurations. Such large amounts of computational time are accept- 
able for proof-of-concept studies and selective analysis. With the advent of the next 
generation of parallel supercomputers, airplane design and optimization using nonlin- 
ear CFD should become routine. An essential element in design and optimization is 
acquiring the sensitivity of functions of CFD solutions with respect to control param- 
eters. For aerodynamic surfaces, the control parameters specify the surface shapes. 
This affects the surface grid and the field grid which, in turn, affects the flow field 
solution. 

Sensitivity analysis (SA) provides a natural systematic means for both an- 
alyzing and predicting the behavior of physical approximations and computational 
systems or for identifying significant input parameters in a system. The system out- 
puts are assumed to be functionally dependent upon the system inputs. The output 
changes in response to specified changes in the input; however, everything within the 
system is normally assumed to be fixed. Changes in the system outputs are related 
to the changes in the inputs through a sensitivity derivative (SD) matrix, or system 
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jacobian. The SD matrix may be used to control processes or designs that depend 
upon the system output. 

Procedures for MDO of engineering systems have been addressed by Sobieski 
and others [97], Sobieski proposes a unified system SA guided by system derivatives. 
Aerodynamics plays a central role which is connected to other disciplines. The ob- 
jective and constraints are provided by the output functions of these several other 
disciplines. Each single discipline is then to supply not only the output functions 
for the constrained optimization process, but also the derivatives of all these output 
functions with respect to its input variables. 

Numerous research efforts have examined the issue of efficient computation 
of SD for CFD. Typical techniques for computing sensitivities include by hand, by 
use of a symbolic expression differentiator and by approximation via divided differ- 
ences. Unfortunately, none of these techniques can be used to deliver fast and reliable 
derivatives in a flexible and timely fashion for large computer codes. Hand coding 
of derivatives is impractical and symbolic approaches may require as much effort as 
hand coding. Divided differences may not be accurate and are obtained too slowly. 

Automatic differentiation (AD) promises to address the need for a flexible 
and scalable technology capable of computing derivatives of large codes accurately, 
irrespective of the complexity of the model. In this study AD has been successfuly 
applied to obtain sensitivity derivatives required for an MDO procedure. Incremental 
iterative form, also known as the “delta” or “correction” form, is another successful 
approach for obtaining the solution state vector from the nonlinear governing flow. 
References [98], discuss the benefits of using this form to solve the large systems of 
linear equations needed to obtain SD. This incremental iterative formulation is very 
flexible and Korivi et al. [99] have demonstrated the use of this strategy to efficiently 
and accurately calculate quasianalytical sensitivity derivatives for a space-marching 
3D Euler code with supersonic flow over a blended wing-body configuration. 
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6.2 Automatic Differentiation and ADIFOR 

Automatic Differentiation (AD) is a collection of computer science tech- 
niques which permit one to automatically calculate the derivatives of information 
generated by a computer program with respect to any parameter intervening in its 
calculation. AD is essentially an automatic implementation of the chain rule of diffren- 
tiation based on tracking the relationships between dependent and independent vari- 
ables. Typically, to calculate the derivative of the output of a program with respect 
to its input, one modifies the original program by insertion of specialized instruc- 
tions which identify relevant independent and dependent variables. The program is 
then modified automatically by a preprocessor which enhances it to calculate deriva- 
tives. The enhanced program is compiled conventionally, linked with special run-time 
libraries (if required) and executed to generate not only the original program’s de- 
pendent variable but also their derivatives with respect to the independent variables. 

There are two modes of AD. In the first, the forward mode, the chain rule 
is evaluated from the input to the output; in this mode, the computational cost in- 
creases with the number of outputs. In this mode, the chain rule is evaluated from 
the output to the input. While it can be much faster than the forward mode, this 
reverse mode can place enormous demands on computer storage and requires special 
memory handling. AD is distinct from finite difference or symbolic manipulation 
techniques. The former, based on perturbations of a program’s input, generates ap- 
proximate derivatives which can be affected by round-off and truncation errors [100]. 
While an exact technique, the later tends to generate very cumbersome expression 
for the derivatives. 

The tool used in the present study is called ADIFOR [101] (automatic differ- 
entiation of Fortran). The tool is jointly developed by Argone National Laboratory 
and Rice University and it differentiates program written in Fortran 77. That is, 
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given a Fortran subroutine (or collection of subroutines) that describe a function and 
an indication of which variables in parameter lists or common blocks corresponds to 
“independent” and “dependent” variables with respect to differentiation, ADIFOR 
produces Fortran 77 code that allows the computation of the derivatives of the depen- 
dent variables with respect to the independent ones. ADIFOR employs a hybrid of 
the forward and reverse modes of AD. That is, for each assignment statement, code is 
generated for computing the partial derivatives of the result with respect to the vari- 
ables on the right-hand side and then the partials are employed in the forward mode 
to propogate overall derivatives. The result is a significant decrease in complexity 
when compared to the forward mode of implementation. The ADIFOR tool produces 
portable Fortran 77 code and accepts almost all of Fortran 77, in particular, arbi- 
trary calling sequences, nested subroutines, common blocks etc. ADIFOR-generated 
code can be used in various ways. Instead of simply producing code to compute the 
Jacobian J, ADIFOR produces code to compute J*S, where the “seed matrix” S is 
initialized by the user. Therefore, if S is the identity, ADIFOR computes the full Jaco- 
bian; whereas if S is just a vector, ADIFOR computes the product of the JACOBIAN 
by a vector. The running time and storage requirements of the ADIFOR-generated 
code are roughly proportional to the number of columns of S, so the computation 
of Jacobian- vector products and compressed Jacobians requires much less time and 
storage than does the generation of the full Jacobian matrix. 


6.3 Theoretical Formulation 


An implicit representation of a physical system can be modeled mathemat- 


ically as 


f{h,g(h)) = o 


( 6 . 1 ) 


where G and H are dependent and independent variables, respectively. The function 
F can have algebraic, differential, integral or integral- differential characteristics. 
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The quantities G and H can be either scalar or vector depending on the nature of the 
physical model. The sensitivity of G with respect to H can be obtained by implicit 
differentiation of Eq. (5.1) 



dF 


dG 



( 6 . 2 ) 


The coefficients, || and [|£], can be obtained, provided that the solution to 
Eq.(5.1) is known. Equation (5.2), now a set of algebraic equations, can be eas- 
ily solved for the sensitivity derivative, If {§77} an< ^ [fg] are not available, a 

finite difference approach can be adopted. The central difference approximation of 
1 can be devised as 


/ dG \ _ G{H + AH) - G(H - AH) 

2 AH 


where AH is a small perturbation of a specified parameter. Although the implemen- 
tation of the finite difference approach is comparatively easy, it has the disadvantage 
of being computationally expensive. Also, the choice of AH is crucial for accuracy 
of the derivative. A large values of AH may lead to inaccurate derivatives while a 
small value may result in round-off errors. 


6.4 Application of ADIFOR to Potential Flow Code 

(PMARC) 

Application of ADIFOR to advanced flow code has been done by Green et 
al.[ll4]. The results have been very encouraging and in the present study ADIFOR 
is applied to the potential flow code PMARC. 

Figure 6.1 indicates an analysis system with input as S and D and output as 
Q. Both S, D and Q may be scalar, vector, or array quantities, and each may involve 
one or more variables. The input S and D for the system consist of several 
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Fig. 6.1 Typical system with ADIFOR applied. 
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types of variables such as those that specify initial and boundary conditions, material 
properties, constraints, physical dimensions and approximations, design variables, and 
numerical solution parameters, etc. Similarly the output Q for the system my consist 
of local and global solution properties, accuracy measures, and performance indicators 
etc. When ADIFOR is applied to the above described system, the output consists of a 
combination of derivatives of the ouput system functions with respect to input system 
variables. These are indicated as dotted line in the figure. Application of ADIFOR to 
Fortran codes requires the specification of the independent and dependent variables 
to be used in forming the SD matrix. 

In this study the computations of CL,CD,andC\i have been used as the 
dependent variables. Application of ADIFOR to PMARC was performed in a very 
simple and straight forward manner. Minor changes to PMARC code was required 
for ADIFOR processing to be accomplished. The PMARC code was passed on to the 
ADIFOR as input. ADIFOR differentiated through the entire solution algorithm, the 
specified dependencies were traced and a new SD code was generated as required. The 
resulting SD modules were then assembled into a working code and the initial results 
were generated quickly. The code was run on an SGI Indigo machine and various 
test cases were examined, and comparisons with direct differentiation procedure have 
been made. 


6.5 Application Of ADIFOR to Grid Generator (RAPID) 

Unlike aerodynamic considerations, the grid sensitivity analysis has been 
used on structural design models for a number of years. In this context, grid sensi- 
tivity can be thought of as perturbation of structural loads, such as displacement or 
natural frequency, with respect to finite element grid point locations [102]. Two basic 
approaches have been cited for grid sensitivity derivatives. The first approach, known 



as implicit differentiation, is based on implicit differentiation of discretized finite el- 
ement system. The other, which is based on the variation of continuum equations, 
is known as variational or material derivative approach. Gradient based techniques 
applied to aerodynamic configurations optimization require the determination of grid 
sensitivity 9 *p r/ )- In the past, in order to evaluate such deriva- 

tives, each expression would have to be differentiated and chain ruled through out 
the mathematical system, either by hand or with the aid of a computer-aided alge- 
braic manipulation system. The simplest way to obtain grid sensitivity is to vary the 
control parameters, one at a time, and finite difference the results. This, however, 
is proven to be computationally inefficient compared to analytical or semi-analytical 
differentiation of the grid equations. Also, the proper choice of a step size is not 
trivial and an improper choice might result in round-off error accumulation. The 
finite difference approach should only be used as the last resort when the extreme 
complexity of the grid equations dictates no other alternatives. For a less compli- 
cated grid equations, a semi-analytical approach would be more appropriate. The 
semi-analytical approach consists of analytical differentiation of the original function 
with respect to an intermediate function, the derivative of which is then evaluated 
numerically. It combines the efficiency of the analytical approach with the ease of 
implementation of the finite difference approach. 

The analytical approach to the grid sensitivity problem is evaluation of the 
grid sensitivity coefficient by direct analytical differentiation of the grid equation. 
For most cases, the grid equation is not directly differentiable, although there are 
schemes that such differentiations are feasible. The algebraic grid generation scheme, 
such as Two-Boundary Grid Generation (TBGG), was successfully differentiated by 
analytical methods by Sadrehaghighi [108] and very accurate results were provided. 
The analytical approach has the advantage of being exact, thus, avoids the round-off 
errors associated with numerical approaches. Due to the time consuming nature and 
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tedious process of analytical approach, in this study grid sensitivity was obtained by 
applying ADIFOR software to the grid generator RAPID. 

Geometric inputs to the RAPID software consists of wing planform and 
wing section definitions, fuselage defining parameters and grid spacing control pa- 
rameters. These parameters are identified as independent variables in ADIFOR, and 
the output surface defining grid coordinates (x,y,z) as the dependent variables. ADI- 
FOR successfully differentiates the entire RAPID software by identifying the various 
dependencies and an enhanced code is generated. The enhanced ADIFOR grid gen- 
erator RAPID 4 £> not only generates the grid, but also generates the derivatives or 
grid sensitivities to the geometric parameters. 


6.6 Optimization Problem 

An objective of a multidisciplinary optimization of a vehicle design is to extremize 
a payoff function combining dependent parameters from several disciplines. Most 
optimization techniques require the sensitivity of the payoff function with respect to 
free parameters of the system. For a fixed grid and solution conditions, the only free 
parameters are the surface design parameters. Therefore, the sensitivity of the payoff 
function with respect to design parameters are needed. 

The optimization problem is based on the method of feasible directions 
[104,105] and the generalized reduced gradient method. This method has the advan- 
tage of progressing rapidly to a near-optimum design with only gradient information 
of the objective and constrained functions required. The problem can be defined 
as finding the vector of design parameters X/j, which will minimize the objective 
function /(Xp) subjected to constraints 


Sb(X D ) <0 j = 1, m 


(6.4) 
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and 

X l D < X D < X U D (6.5) 

where superscripts denote the upper and lower bounds for each design parameter. 
The optimization process proceeds iteratively as 

X n D = X n D ~ l +7S n (6.6) 

where n is the iteration number, S n the vector of search direction, and 7 a scalar 
move parameter. The first step is to determine a feasible search direction S n , and 
then perform a one-dimensional search in this direction to reduce the objective func- 
tion as much as possible, subjected to the constraints. 

The present optimization strategy is based on maximizing the lift coefficient, 
Cl, in response to surface perturbation, subject to pre-determined design constraints. 
Upper and lower bounds set for each design parameter and the sensitivity derivatives 
of the objective function, §§^1 and the constraint, are obtained as previously 

described. Throughout the analysis, the drag coefficient, Cp, is to be no greater than 
the value of the initial design. The strategy, illustrated in Fig. 6.2, requires that the 
grid and grid sensitivity derivatives be provided dynamically during the automated 
optimization process. 












Chapter 7 


RESULTS AND DISCUSSION 

Two test cases are considered to demonstrate the feasibility of current pro- 
cedure. For each case, the grid and flow sensitivity coefficients of the field have been 
obtained. The sensitivities of the total forces (i.e., Lift and Drag coefficients) are 
tabulated for optimization purposes. The first test case, a symmetrical generic air- 
plane with 14 surface defining parameters (Fig. 7.1), has been used mainly to exhibit 
the accuracy of grid sensitivity coefficients with those obtained using finite differ- 
ence approach before proceeding to a realistic configuration. The second test case, 
a HSCT type configuration (Fig. 7.2), has been used to extend the analysis to do 
a three dimensional optimization. An optimization module has been integrated into 
the overall procedure to optimize the geometry using the resultant sensitivity coef- 
ficients. The improved design is used for the Euler study where an Euler type two 
block volume grid is constructed using GRIDGEN software and solutions obtained 
using the TLNS3D code. The CSCMDO[106] software is also used to transform the 
grid from the original geometry to the new optimized geometry. 

7.1 Grid Sensitivity 

Grids obtained from RAPID software, shown in Fig. (3.1), is considered for 
grid sensitivity analysis. Grid sensitivity study was performed on the HSCT type 
configuration shown in Fig. (7.2) with fourteen surface defining design variables. 
The surface grid sensitivity with respect to the vector of design parameters, Xp, is 
obtained from the ADIFOR differentiated code RAPID. 
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Fig. 7.1 Generic airplane for case 1. 
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Figure 7.3 shows the x-coordinate sensitivity with respect to chamber. The 
highest contour levels are, understandably, located at the point of maximum cham- 
ber. Since the wing sections are defined as a NACA four-digit wing sections, this is 
positioned at about 0.3 of the chord length from the leading edge[ 107] . The positive 
and negative contour levels correspond to the upper and lower surfaces. The sensi- 
tivity levels decrease when moving away from the location of maximum chamber. 

Typical CFD calculations are performed on a computational mesh that is 
“body-oriented”. Changes in the geometric shape results in the movement of grid 
points throughout the entire mesh. The benchmark for comparison of these grid sen- 
sitivity terms is by performing finite difference. If forward difference approximations 
are selected, for example, the mesh generation code is used to produce one additional 
perturbed grid for a slightly perturbed value of geometric shape design variable which 
in this case is the camber. Finite difference at each of the grid coordinate is calculated 
and the result is shown in Fig. 7.4. It is seen that the results obtained by AD IFOR 
and the finite difference are in very good agreement, thus confirming the accurate 
results obtained by ADIFOR. 

Having confirmed the results obtained from ADIFOR, and to further evalu- 
ate the results, next the design variable chord is choosen as the independent param- 
eter. Figure 7.5 shows the x-coordinate sensitivity with respect to the chord. The 
concentration of contours are near the leading edge, in the x-direction of the wing 
fuselage intersection thus confirming the maximum effect at that point. The increase 
in chord moves the tip of the wing fuselage airfoil section towards the nose of the 
airplane configuration while the trailing edge is kept fixed. Figure 7.6 



Fig. 7.3 X-coordinate sensitivity with respect to camber. 



Fig. 7.4 Finite difference X-coordinate sensitivity with 
respect to camber. 



Fig. 7.5 X-coordinate sensitivity with respect to chord 



Fig. 7.6 Y-coordinate sensitivity with respect to chord. 



Fig. 7.7 Z-coordinate sensitivity with respect to chord. 
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shows the y-coordinate sensitivity with respect to chord. Unlike the x coordinate 
sensitivity the contour levels this time is concentrated in the y-direction indicating a 
maximum change. Figure 7.7 shows the z coordinate sensitivity with respect to chord. 


7.2 Flow Sensitivity 

When AD is applied directly to the potential flow code PMARC, the re- 
sulting AD-enhanced code calculates the required sensitivity derivatives through an 
iterative process. The ADIFOR procedure generates a new version of the potential 
flow code that has the capability to calculate the derivatives of lift, drag, and pitching 
moment with respect to a wide variety of different types of input parameters (includ- 
ing parameters related to the geometric design). 

Both geometric and non-geometric design variables are considered to evalu- 
ate the accuracy of ADIFOR enhanced PMARC (P M ARC ad)- Angle of attack (a) 
is considered as the non-geometric design variable and both sensitivities of lift and 
drag are computed. For the geometric design variable, wing thickness and fuselage 
diameter are considered. The values are compared with the finite difference results 
and are tabulated in Table 7.1. It is seen that the values obtained by ADIFOR are 
in good agreement with the finite difference values, thus confirming the successful 
differentiation of the PMARC code. 

7.3 Optimization Problem 

An objective of multidisciplinary optimization of a vehicle design is to ex- 
tremize a payoff function combining dependent parameters from several disciplines. 
Most optimization techniques require the sensitivity of the payoff function with re- 
spect to free parameters of the system. For a fixed grid and solution conditions, 
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Table 7.1. Comparison Of ADIFOR results with finite difference 
for geometric and nongeometric design variables. 


|§| 

ADIFOR 

Finite Diff. 


$C* , 

.. 




| Angle Of Attack 

lllilil 

lllKiil 

019974 

-*01127 

WingTMdoutsi 

0171446 

’ : ;• ‘ •' : . : x ; - . • ij -• 

-Qxizm | 

0171#; 

. | :£j:§ 

|gX07547 

Fuselage Blanwter 

Siiftp 

-0.09421 

' ~219fla|i: 

ffao#ai 



LOT 


the only free parameters are the surface design parameters. Therefore, the sensitivity 
of the payoff function with respect to design parameters are needed. 

The present optimization strategy is based on maximizing the lift coefficient, 
Cl, in response to surface perturbation, subject to pre-determined design constraints. 
Upper and lower bounds are set for each design parameter and the sensitivity deriva- 
tives of the objective function, fj£-, and the constraint, are obtained as pre- 

viously described. Throughout the analysis, the drag coefficient, Co , is to be no 
greater than the value of the initial design. The strategy, illustrated in Fig. 6.2, 
requires that the grid and grid sensitivity derivatives be provided dynamically during 
the automated optimization process. 

Optimization of the HSCT type configuration shown in Fig. 7.2 was carried 
out on SGI machine with memory capacity of 512 MB. Sixteen design variables were 
selected for the optimization process. A total of twelve design optimization cycles 
were performed and each iteration took approximately 7.5 min of cpu time. It was 
noted that the lift which was initially 0.01712 became 0.0748. The initial and final 
shapes with shaded Cp plots are shown in Figs. 7.8 and 7.9. The comparison of 
the two shapes, before and after the optimization cycle is shown in Fig. 7.10. A 
considerable increase in the length of the inner and outer wings with an increase in 
the wing planform area is seen. A decrease in the wing and fuselage thickness with 
an increase in camber is also noted. 


7.4 Euler Flow Solutions 

To verify the results obtained from the optimizer, it was suggested to perform 
Euler calculations over the initial and final shapes of the HSCT configuration. A 
semidiscrete, cell-centered finite volume algorithm TLNS3D , based on a Runge-Kutta 
time-stepping scheme, is used to obtain the Euler solutions around the HSCT 



Fig. 7.8 Cp over the initial HSCT configuration. 



Fig. 7.9 Cp over the final HSCT configuration. 
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Fig. 7.10 Comparison of initial and final shapes. 
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type configuration. The efficiency of the numerical scheme is greatly enhanced by 
taking advantage of the multigrid acceleration technique. A two block, C-0 mesh 
with 142x82x42 grid, shown in Fig. 4.12, is used to obtain converged solution at 
Mach number of 2.4. Figure 7.11 shows the Cp plot over the original geometry. A 
lift of 0.01748 was obtained, and it compared very well with the potential flow case. 

Next a volume grid was generated over the optimized HSCT surface shown 
in Fig. 7.9. In this case the CSCMDO software was used to interpolate points from 
the volume grid generated over the original geometry, thus saving time by not gener- 
ating the grid from scratch. Figure 7.12 shows the shaded Cp plot over the optimized 
configuration. Figures 7.13 and 7.14 show the line plot of Cp at the crank of the 
wing. It is noted that the distribution of pressure is very smooth and well behaved 
over the surface of the optimized configuration, thus confirming the trend of the re- 
sults obtained from the optimizer. 



Fig. 7.12 Euler flow on the optimized configuration 
derived from potential flow optimization. 






Chapter 8 


CONCLUSIONS AND RECOMMENDATIONS 

An algorithm is developed to define the surfaces of aerodynamic configura- 
tions for design optimization. The two schemes investigated are Non-Uniform Ratio- 
nal B-Splines (NURBS) and Partial Differential Equation (PDE). NURBS parametriza- 
tion defines the surface by a set of ordered control points. These control points act 
as a set of design parameters which is used in an optimization process. The PDE 
technique offers a unique type of parametrization where a surface is defined by a 
fourth order partial differential equation. This procedure generates a blend surface 
between two sets of curves. The design parameters in this case is the constants used 
in the definition of the two curves. 

The PDE technique is used towards the parametrization of a HSCT type 
configuration. Inclusive in this definition are surface grids, volume grids, and grid 
sensitivity. The design variables are incorporated into the boundary conditions, and 
the solution is expressed as a Fourier series. The fuselage has circular cross section, 
and the radius is an algebraic function of four design parameters and an independent 
computational variable. Volume grids are obtained through an application of the 
Control Point Form method. 

A graphic interface software is developed to represent the PDE surface. The 
software has the capability to dynamically change the surface with the change in input 
design variables. Various options and features are provided to enhance the quality of 
the software which gives a competitive outlook. 

Grid sensitivity with respect to surface design parameters and aerodynamic 
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sensitivity coefficients based on potential flow is obtained using an Automatic Dif- 
ferentiation precompiler software tool ADIFOR. Aerodynamic shape optimization of 
the complete aircraft with twenty four design variables is performed. Unstructured 
and structured volume grids and Euler solutions are obtained to demonstrate the 
feasibility of the new surface definition. 

Future investigations should include the implementation of present approach 
using larger grid dimensions, adequate to resolve full physics of viscous flow analysis. 
A grid optimization mechanism based on grid sensitivity coefficients with respect to 
grid parameters should be included in the overall optimization process. An optimized 
grid applied to present geometry, should increase the quality and convergence rate of 
flow analysis within optimization cycles. Other directions could be establishing a link 
between the graphic software and the optimization procedure, such that the changes 
in the design variables in each optimization cycle is visible to the user. Another 
contribution would be the extension of the current algorithm to represent complex 
configurations. A hybrid approach can be selected where certain sections or skeletal 
parts of a surface are specified analytically and interpolation formulas are used for 
intermediate surfaces. 



